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ABSTRACT 


An  analysis  of  the  optimum  statistical  restoration  of 
quantized  signals  is  presented.  The  restoration  is  based 
upon  minimizing  the  mean  square  error  between  the  input  to 
a quantizer  and  its  estimate.  Since  a quantizer  is  a 
nonlinear  device,  the  estimation  equation  which  is  derived 
achieves  an  optimum  nonlinear  restoration.  Its  solution 
requires  complete  statistical  knowledge  of  the  quantizer 
input.  Available  statistical  information  usually  includes 
the  marginal  distribution  of  each  of  the  input  variables 
and  the  correlation  between  them.  Hence  a technique  is 
developed  for  generating  correlated  multidimensional 
probability  density  functions  based  on  this  available 
information.  The  technique  is  applied  to  gaussian, 
laplacian,  and  Rayleigh  density  functions.  These 
multidimensional  density  functions  characterize  the  outputs 
of  transform  coders,  DPCM  coders,  and  PCM  coders, 
respectively.  The  quantized  outputs  of  these  coders  are 
then  restored  by  utilizing  the  multidimensional  densities 
in  the  estimation  equation.  Examples  of  images  which  have 
been  coded  and  restored  by  these  techniaues  are  presented. 
The  results  reveal  a mean  square  error  reduction.  To 
achieve  a visually  subjective  improvement  also,  a weiqhted 
mean  square  error  criterion  is  employed,  where  the 
weighting  corresponds  to  characteristics  of  the  human 
visual  system. 
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CHAPTER  1 
INTRODUCTION 


Quantization  is  a process  inherent  in  all  digital 
systems.  Basically,  quantization  occurs  whenever 
continuous  physical  properties  are  represented  numerically. 
When  this  representation  takes  place  in  a digital  computer, 
the  quantization  effect  is  called  round-off,  ot  truncation. 
As  a mathematical  operation,  quantization  is  the  processing 
of  continuous  functions  to  give  a stepwise  output,  or  the 
processing  of  sampled  functions  to  give  a sampled  output. 
Even  the  value  obtained  in  measuring  a continuous  quantity 
is  the  consequence  of  quantization.  But  wherever  the 
occurrence,  a fundamental  aspect  of  quantization  is  that  it 
results  in  an  indeterminacy  and  a lack  of  complete 
information  about  the  particular  property  under 
consideration . 

In  this  report  quantization  is  assumed  to  be  a 
nonlinear  operation  which  occurs  within  a cuantizer — a 
zero-memory  device  that  assigns  an  input  to  one  of  a 
countable  number  of  possible  output  regions.  This  defines 
a broad  class  of  devices  that  includes  coders,  digital 
transducers,  and  analog-to-digi tal  converters.  For  many  of 
these  quantizing  devices  the  input  is  a continuous 
variable;  restricting  a continuous  input  to  a particular 


t 
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region  destroys  some  of  the  information  about  that  input. 
For  discrete  inputs,  combining  their  input  regions  into 
larger,  and  hence  fewer  allowable,  output  regions  also 
decreases  the  amount  of  information  available.  These  are 
both  irreversible  operations  and  the  lost  information 
cannot  be  recovered.  It  is  thus  important  to  optimize  the 
quantizing  process  so  that  this  lost  information  is  held  to 
a minimum. 

The  lost  information  can  be  minimized  by  decreasin', 
the  size  of  the  output  regions  while  simultaneously 
increasing  their  number.  Unfortunately,  this  is  not  always 
possible  or  practical.  It  would  also  oppose  one  of  the 
benefits  of  quantization:  a smaller  number  of  output 
regions  requires  less  processing  and  less  storage.  A 
balance  thus  must  be  attained  between  accuracy  and  economy. 
This  balance  can  be  determined  by  an  analysis  of  the 
quantization  nrocess. 

Most  analyses  of  auantization  to  date  have  focused  on 
just  one  aspect  of  the  problem,  i.e.,  finding  the  best 
quantizing  device  to  minimize  the  information  lost, 
however,  the  eventual  use  of  any  quantized  output  is  to 
accur  'ly  represent  a continuous  signal  input.  The  output 
regions  are  ultimately  utilized  to  estimate  and  restore  the 
original  quantizer  input.  It  is  this  quantization 
restoration  problem  that  has  heretofore  been  neglected  and 
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is  the  subject  of  this  dissertation. 


The  simolest  restoration  procedure  is  to  choose  the 
midpoint  of  each  quantization  interval  as  the  estimated 
value  of  the  original  input.  However,  this  estimate  can  be 
improved,  since  it  is  based  only  on  the  outDut  reqions  of 
the  quantizer.  The  restoration  to  be  described  herein  is 
based  also  on  a priori  knowledge  of  the  statistics  of  the 
quantizer  input.  The  input  is  assumed  to  consist  of 
samples  from  a continuous  random  process.  (The  samplinq 
presents  few,  if  any,  restrictions  because  digital  systems 
reauire  sampled  and  quantized  signals,  and  the  operations 
of  sampling  and  quantizing  are  commutable.)  The  necessary 
statistics  are  the  amplitude  probability  density  function 
of  the  input  samples.  Where  the  complete  statistics  are 
unknown,  a functional  form  for  them  is  developed  from  known 
correlation  functions  and  one-dimensional  distributions. 
This  statistical  information  is  then  combined  with 
knowledge  of  the  quantizer  output  to  provide  an  optimum 
restoration.  The  restoration  is  optimum  with  respect  to  a 
desired  error  criterion. 


A second  problem  requiring  a similar  solution  occurs 
when  quantized  signals  are  transmitted  through  a noisy 
channel.  Because  of  the  errors  that  accrue  during 
transmission,  what  is  received  does  not  exactly  correspond 


to  the  quantizer  output  regions.  The  channel  output 


3 


instead  equals  the  quantizer  output  only  within  a specified 
probability.  To  achieve  an  optimum  restoration  the 

following  available  information  must  be  utilized: 

1.  the  a priori  distribution  of  the  quantizer  input 

2.  the  structure  of  the  quantizer 

3.  the  transmitted  quantizer  output  reqion  (which  may 
be  in  error) 

4.  the  channel  error  structure 

The  existence  of  the  last  constituent  induces  a 
modification  in  the  solution  discussed  previously.  The 
modification  is  also  considered  in  this  dissertation. 


In  this  dissertation,  solutions  are  presented  for  the 
restoration  of  quantized  samples  based  on  a priori 
knowledge  of  the  multivariate  probability  density  function 
of  the  quantizer  input.  The  two  cases  considered  are: 
(1) , the  quantizer  output  region  is  known  exactly;  and  (2), 
the  quantizer  output  is  transmitted  over  a noisy  channel 
and  hence  not  known  exactly. 


These  two  situations  arise  in  the  coding  and 
transmission  of  images.  Quantizers  are  an  integral  part  of 
all  image  coding  systems.  The  goal  of  these  systems  is  to 
make  a coded  image  as  similar  as  possible  to  an  original 
image.  Unfortunately  cost,  complexity,  and  hardware 
constraints  often  force  a suboptimal  coding  scheme  which 
results  in  a degraded  image.  The  application  of  the 
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quantization  restoration  techniques  discussed  above  can 
improve  these  degraded  images.  Experimental  verification 
of  this  improvement  is  obtained  by  restoring  images  which 
have  been  coded  and  quantized. 


CHAPTER  2 


HISTORY  OF  QUANTIZATION  AND  SIGNAL  RECONSTRUCTION 


Developments  in  quantization  have  closely  paralleled 
advances  in  digital  systems.  Although  research  had  been 
:onducted  into  areas  such  as  uniform  statistical  grouping 
(which  may  be  considered  quantization)  as  early  as  1898 
1 1 ] , it  is  only  since  1947  that  deliberate  attempts  have 
been  made  to  understand  the  process  of  Quantization.  By 
1947  vacuum-tube  technology  had  reached  the  stage  for  which 


digital  systems 

were  both  possible  and  practic 

al.  At 

that 

time  the 

concepts  and  the  value 

of  PCM 

(pulse 

code 

modulation 

— the 

first  major  application 

of 

both 

aunatization  and 

digital  hardware) 

were 

just 

being 

discovered 

and 

made  known  [2,3,4] 

. Bennett  [5] 

then 

undertook 

an  intensive  investigation 

of  the 

spectra  of 

quantized  signals.  Bennett  analyzed  uniform  quantizers, 
such  as  those  uitlized  in  PCM  systems,  and  found  the 
characteristics  of  their  output  spectra  for  a white  noise 
input  spectrum.  Since  Bennett's  initial  work,  developments 
in  quantization  have  proceeded  along  three  basic  lines: 

1.  analyzing  the  results  and  the  process  of 

quantization 

2.  optimizing  the  quantization  process  with  respect 
to  various  criteria  and  goals 

reconstructing  quantized  signals  to  minimize  the 
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3. 


deq r adation  incurred  throuqh  ouantization 
Each  of  these  cateqories  will  be  considered  separately  in 
the  followinq  paragraphs,  and  siqnificant  developments  will 
be  discussed  in  rouqhly  chronoloq ical  order. 

2.1  Anal /sis  of  Quantization 

One  of  Bennett  s conclusions  was  that  quantization 
uncertainty  or  noise,  for  a quantizer  with  many  levels , i s 
uniformly  distributed  throuqhout  the  signal  band.  This 
result  was  supported  by  Widrow  in  1956  [6]  in  studies  of 


h. 

k 

P',, 


the  probability  density  functions  of  quantized  signals  and 
quantization  noise.  Widrow  concluded  that  a quantizer 
could  be  modelled  as  a source  of  uniform,  independent 
noise.  In  a later  paper  [7],  Widrow  attempted  to  define 
the  limits  of  the  region  over  which  his  additive  noise 
conclusion  would  be  valid,  and  then  extended  the 
statistical  results  to  two  dimensions  (i.e.,  the 
quantization  of  two  correlated  samples) . Myers  [8] 
extended  Bennett's  analysis  to  the  case  of  the  uniform 
quantization  of  a signal  corrupted  by  gaussian  noise  and 
derived  the  resultant  noise  distribution,  which  is  no 
longer  uniform.  Velichkin  [9] , calculated  the  correlation 
function  and  output  spectrum  of  a quantizer.  In  addition, 
Velichkin  considered  the  more  general  cases  of  nonuniform 
quantization  levels  and  arbitrary  input  signal  spectrums. 


Velichkin 's  results,  for  the  gaussian  case,  were  in  the 


in:. 


— * — — 


form  of  an  infinite  summation  of  Hermite  polynomials  of 
increasing  order  which,  unfortunately,  cannot  be  evaluated 
without  simplifying  assumptions. 


Robertson  [10]  surmounted  this  difficulty  by 
evaluating  combinations  of  the  terms  of  the  summation  such 
that  the  combinations  tend  to  zero  and  hence  obtained 
output  spectrums  for  nonlinear  and  nonsymmetrical 
juantizers  and  for  arbitrary  input  spectrum  shapes.  Chan 
and  Donaldson  [11]  obtained  a further  generalisation  by 
finding  the  correlation  function  and  spectrum  of  a 
quantized  gaussian  signal  transmitted  over  a discrete 
memoryless  channel.  Their  results  reduce  to  those  obtained 


by  Velichkin  when  the  channel  is  noiseless.  For  very 
coarse  quantization,  Curry  and  Vander  Velde  [12]  suggested 
modelling  the  quantizer  as  a gain  element,  whose  value  is 
equal  to  the  random-input  describing  function,  plus  an 
additive  noise  source.  The  inclusion  of  a gain  element 
causes  the  noise  source  to  appear  more  nearly  white.  The 
quantization  process  can  then  be  analyzed  more  easily. 

2.2  Optimum  Quantizers 

The  analyses  of  q -.ntization  described  above  were  all 


based  on  the  concept  that  quantization  introduces  a noise 
or  distortion.  A number  of  researchers  have  attacked  the 

fun  H 1 nrnhl  oin  f minimi  a m c*  ^ i f 


distortion  by  varying  the  location  of  N quantizing  levels, 
given  the  characteristics  of  the  quantizer  input.  They 
attempted  to  do  this  in  an  optimal  fashion  under  different 
assumptions  and  conditions,  and  according  to  various 
cr iter ia . 

In  1951  Panter  and  Dite  [13]  tried  to  minimize  the 
mean-square  quantization  error  by  utilizing  statistical 
properties  of  the  signal.  They  developed  an  optimum 
nonuniform  quantizer  based  on  the  following  assumptions: 

1.  the  quantizer  is  symmetrical  about  zero 

2.  the  probability  density  function,  p(x) , is  an  even 
function  and  is  constant  over  each  quantization 
interval 

3.  the  signal  is  limited  to  the  range  |-V,V] 

4.  a signal  quantized  to  a particular  int  rval  is 
restored  to  the  midpoint  of  that  interval 

The  resultant  quantizer  restoration  levels,  rk , can  then  be 
calculated  from 


■/ 


2kV/N  -1/3 

[P ( x ) ] dx 

0 


V 


j>(x 


-1/3 

)]  dx 


(2.1) 


for  k=l , 2 , . . . , n and  where  N=2n+1  is  the  total  number  of 
levels.  The  total  distortion  power  for  this  choice  of 


j f±~l. 


9 


levels  is 


£ = — o t / t P ( x ) 1 dx)  (2.2) 

3N  J 0 

Panter  and  Dite  also  suggested  that  this  nonuniform  level 
spacing  could  be  realized  by  "companding" — compressing  the 
original  signal  by  a nonlinear  function  such  as  a 
logarithm,  performing  a uniform  quantization,  and  then 
expanding  the  result  by  means  of  the  inverse  of  the 
nonlinear  function.  This  type  of  system  was  later  analyzed 
by  Smith  [14]  who  provided  a method  for  choosing  the 
parameters  of  the  nonlinearity  with  respect  to  a 

mean-square  error  criterion. 

In  a fundamental  paper  in  1960,  Max  [15]  derived  the 
necessary  equations  for  finding  the  parameters  of  a 
quantizer  having  minimum  distortion  with  respect  to  a 

convex  error  criterion.  For  a fixed  number  of  quantization 
intervals,  N,  the  decision  levels,  d^ , and  the  restoration 
levels,  r^ , are  obtained  by  a recursive  solution  of 

dk=(rk  + rk-l)/2  k = 2 , 3 , . . . , N (2.3a) 


[\+i 

J \ 


be  (x-r  ) 

P ( x ) dx=0 

b* 


k=l , 2 , . . . , N 


(2.3b) 
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where  p(x)  is  the  distribution  function  of  the 
input  and  e(*)  is  a convex  error  function, 
these  equations  do  not  require  o(x)  to  be  constan 
range  of  the  quantizer.)  In  this  notation, 
quantized  to  the  interval  (c]^  , c|^  ) would  be  re 

the  point  r . The  set  of  simultaneous  equations 
in  eq.  2.3,  except  for  trivial  cases,  cannot 
explicitly  and  so  must  be  evaluated  numerically, 
derived  the  equation  which  provides  the  optimum 
parameters  for  a uniform  spacing  of  levels,  and 
both  the  optimum  uniform  and  optimum 
quantization  levels  for  a gaussian  orobabili 
function  and  a mean-square  error  criterion. 


quanti zer 
(Note  that 
t over  the 
a signal 
stored  to 
contained 
be  solved 
Max  also 
quantizer 
tabulated 
nonuni  form 
ty  density 


The  difficulty  in  solving  eq.  2.3  explicitly  has  led  ? 
number  of  researchers  to  consider  various  approximations. 
Garmash  [16]  simplified  the  integrals  in  eq . 2.3b  by  the 
trapezoidal  rule  and  reduced  the  number  of  equations  by 
choosing  the  size  of  the  smallest  interval.  His  results 
are  valid  only  for  a finite  signal  range.  Roe  [17] 
approximated  the  probability  density  function,  p(x),  by  the 
first  two  terms  of  its  Taylor  series  expansion  about  the 
midpoint  of  each  quatization  interval  to  obtain 


l/(t+l 

[p(x) ] 


dx  « 2ak+b 


(2.4) 
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toe  k=l ,”,... ,N-1  where  t is  the  order  of  the  error 
criterion  ( t=  2 for  a menn-sauare  error)  and  a and  b are 
constants.  This  relation  approximately  provides  the 
decision  levels  (and  restoration  levels  for  k'  = k+l/2)  for 
differentiable  probability  density  functions,  but  the 
resultant  quartizer  is  identical  to  that  obtained  by  Panter 
and  Dite  [13].  Algazi  [18],  ir.  attempting  to  find  some 
simpler  suboptimal  algorithms,  also  rederived  eq.  2.2. 
Williams  [19]  published  a closed-form  solution  to  eq.  2.3 
for  the  special  case  of  a laplacian  distribution,  but  erred 
in  assuming  the  restoration  point  to  be  the  median,  rather 
than  the  centroid,  of  the  quantization  interval.  The 
correct  quantization  levels  for  a laplacian  (and  also  a 
gamma)  distribution  were  later  calculated  numerically  [20]. 


No  matter  which  method  is  utilized  to  calculate  the 
quantization  parameters,  however,  it  will  fail  if  the 
probability  distribution  is  such  that  a unique  minimum 
distortion  point  does  not  exist.  For  a mean-square 

distortion  measure,  Fleischer  [21]  derived  the  sufficient 
conditions  under  which  a unique  optimum  quantizer  can  be 
found.  To  this  end  also,  Bruce  [22]  used  dynamic 

programming  to  find  optimum  quantizers  for  a variety  of 
convex  error  criteria,  and  checked  their  uniqueness  by 
locating  and  comparing  all  other  extrema  of  the  error 
surfaces . 


Since  quantizers  usually  operate  on  sampler]  siqnals, 
Velichkin  (23)  and  Goodman  [24]  considered  the  joint 


optimization  of  both  sampling  and  quantizing.  Velichkin 


calculated  the  optimum  parameters  to  achieve  a minimum 


mean-square  error,  and  Goodman  compared  this  result  to  a 


lower  bound  obtained  from  ra'-e  distortion  theory. 


It  has  been  found  [25]  that  uniform  quantizers 


approach  this  lower  bound,  i.e.,  uniform  quantizers 


asymptotically  have  the  lowest  output  entropy.  Wood  [26] 


and  O'Neal  [27]  have  taken  advantage  of  this  fact  to  derive 


(approximately)  minimum  entropy  quantizers.  Wood  has  shown 


that  for  a fixed  output  entropy,  uniform  quantizers  have 


lower  mean-square  error  than  nonuniform  (Max)  quantizers. 


However,  to  achieve  this  error  reduction,  the  quantizer 


output  must  be  optimally  coded  (i.e.,  with  a Huffman 


variable-length  code)  and  this  causes  buffering  problems. 


This  difficulty  may  be  partially  overcome  by  permutation 


encoding  [28],  but  this  coding  technique  has  the  limitation 


of  requiring  very  long  block-lengths.  For  a non-buffered 


coding  scheme,  a quantizer  having  maximum  output  entropy 


could  be  considered  to  be  optimum.  It  has  been  shown  that 


the  quantizers  with  minimum  average  error  (such  as  Max's) 


are  the  same  as  those  with  maximum  entropy,  within  a 


multiplicative  constant  [29]. 


Maximum  entropy,  or  minimum  error,  quantizers  must  be 


modified  when  their  output  is  transmitted  over  a noisy 


channel  [30].  For  a mean-square  error  criterion,  the 
decision  and  reconstruction  levels  for  a noisy-channel 
quantizer  can  be  found  from 


k=2,3,...,N  (2.5a) 


N r d.  . 

E / xp(x)dx 

i = l J d 

r k-1 , 2 , . . . , N (2.5b) 

N rd.±1 

Z Pik  / 1 1P(x)dx 

i=l  ^ J d. 


where  P,  . is  an  element  of  the  channel  matrix,  P.  For  a 
noiseless  channel  (P=J[)  , eq.  2.5  reduces  to  eq.  2.3. 
Minimum  error  quantizers  also  must  be  modified  for 
nonstationary  inputs  or  for  correlated  input  samples. 
Golding  and  Schultheiss  [31]  and  Stroh  and  Boorstyn  [32] 
presented  ad  hoc  adaptive  quantizers  designed  to  handle 
this  situation.  An  optimum  adaptive  quantizer  has  not  yet 
been  successfully  derived,  nor  has  an  optimum  quantizer  for 
correlated  signals. 
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2.3 


Quantized  Siqnal  Reconstruction 


An  important  area  of  research  in  recert  years  has  been 
the  reconstruction  of  quantized  si7~als.  Some  researchers 
have  attempted  to  restore  the  input  signal  samples  which 
were  degraded  by  a quantizer,  while  others,  realizing  that 
many  quantizers  operate  on  sampled  versions  of  analog 
signals,  have  tried  to  restore  directly  the  analog 
waveform.  In  this  latter  category  are  the  techniques  of 
Ruchkin  [33],  Katzenelson  [341,  Steiglitz  [35],  Goblick 
[36],  Kellogg  [37],  Hayes  [38],  and  Chan  and  Donaldson 
[39].  Each  of  their  approaches  differed  in  the  assumptions 
they  made  ibout  the  effect  of  the  quantizer.  Ruchkin  found 
the  best  (with  respect  to  mean-square  error)  linear  filter 
to  restore  a quantized  and  sampled  gaussian  signal,  under 
the  assumption  that  a quantizer  adds  white  gaussian  noise. 
Katzenelson  also  found  the  best  linear  filter,  but  assumed 
the  signal  is  a sample  from  a Markov  process  and  the 
quantizer  noise  is  colored.  The  same  assumption  about  the 
quantizer  was  made  by  Steiglitz  who  analyzed  specific, 
ncnoptimum,  reconstruction  filters  for  special  input  power 
spectra  and  compared  the  results  to  a rate  distortion 
bound.  Steiglitz  provided  a trade-off  betveen  sampling 
rate  and  number  of  uniform  quantization  levels  for  the  case 
of  a minimum  mean-square  error  and  fixed  capacity. 
Goblick,  on  the  other  hand,  considered  entropy  coding  and 
provided  a trade-off  between  mean-square  error  and  output 


entropy.  A prefilter  was  derived  by  Kellogg  for  the  same 
case  that  Ruchkin  considered,  under  the  premise  that  a 
designer  can  tailor  input  and  output  signals  to  a fixed 
quantizer  by  means  of  prefii.ters  and  postfilters.  Kellogg 
used  an  optimum,  nonuniform  quantizer  for  his  simulations, 
and  presented  a numerical  solution  for  the  linear  filters. 
Hayes  added  a constraint  on  the  signal-to-noise  ratio  for 
his  solutions,  but  otherwise  numerically  analyzed  the  same 
case  as  Kellogg.  An  exact  analytical  method  for  jointly 
optimizing  the  prefilter,  postfilter,  quantizer,  and 
sampling  rate  was  presented  by  Chan  and  Donaldson.  Their 
optimum  prefilter  "whitens"  the  signal  and  removes  the 
redundancy  from  the  resultant  signal  samples.  After 


quantization  and  transmission  occur,  a postfilter  restores 
the  continuous  signal  and  removes  the  distortion  added  by 
both  the  quantizer  and  the  channel.  This  method  was 
general  in  that  it  considered  the  correlation  between  the 
signal  and  noise,  included  the  effect  of  channel  errors, 
and  placed  no  restrictions  on  either  the  signal  and  noise 
spectra  or  the  filter  passbands  <*>. 

Curry  [41]  considered  the  problem  of  finding  linear 
estimates  based  on  quantized  samples.  He  showed  that  these 


<*>  It  should  be  mentioned  that  an  excellent  qeneral 
treatment  of  quantizers,  receivers,  and  noisy  chunnels  is 
available  in  Fine  [40]. 


estimates  can  be  computed  in  two  steps:  (1)  find  the 
conditional  mean  of  the  quantizer  output,  and  (2)  pass  this 
mean  throuqh  the  linear  filter  that  would  have  been  used 
had  the  quantizer  not  existed.  For  a minimum  mean-square 
error  criterion  this  linear  filter  was  found  to  be  the 
Kalman  or  Wiener  filter.  The  conditional  mean  for  a 
gaussian  process  was  found  to  be  approximately 

r=(I-DC_1) (b+a)/2  (2.6) 

where  I is  the  identity  matrix,  C is  the  covariance  matrix 
of  the  quantizer  input  x,  x is  quantized  to  the  interval 
[a,b]  bv  an  arbitrary  quantizer,  and  D is  the  diagonal 
matrix  having  elements 


D= 


2 


(b.  -a.  ) 
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*2.7) 


Eauation  2.6  provides  a minimum  mean-saua  e error  estimate 
of  a quantizer's  input,  based  only  on  its  output,  that  is 
valid  whenever  the  quantization  is  very  fine.  Clements  and 
Haddad  [42]  derived  a recursive  techinque  for  finding  this 
same  result  that  is  also  applicable  to  the  problem  of 
nonlinear  estimation  using  quantized  data. 
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CHAPTER  3 

NONLINEAR  ESTIMATION  WITH  QUANTIZED  MEASUREMENTS 

As  the  previous  chapter  has  shown,  most  analyses  of 
quantization  have  focused  on  the  quantizing  process  itself 
and  on  the  determination  of  an  optimum  quantizer  for  each 
input  sample.  Unfortunately,  little  effort  has  been 
expended  towards  undoing  the  effects  of  quantization  and 
recovering  the  original  input.  The  core  of  the  problem  is 
that  quantizers  treat  each  signal  individually  and  neglect 
the  random  process  which  models  the  signals.  As  the 
analysis  below  shows,  an  optimum  restoration  of  auantized 
signals  must  utilize  the  joint  probability  density  of  the 
original  input  signals. 

3.1  Vector  Quantization 

The  input  to  a quantizer  may  be  either  a continuous 
function  or  discrete  samples.  A quantizer  processes  a 
continuous  function  to  give  a stepwise  continuous  output  or 
processes  a sampled  function  to  give  a sampled  output.  All 
digital  systems  require  signals  which  are  both  sampled  and 
quantized.  Sampling  and  quantizing  are  mathematically 
commutable  operations,  i.e.,  the  result  is  the  same  whether 
a signal  is  first  quantized  and  the  resultant  stepfunction 
is  sampled,  or  if  the  signal  is  first  sampled  and  then  the 
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samples  are  quantized.  Hence  in  this  analysis  the  input  to 
a quantizer  is  assumed  to  be  a sampled  random  process. 

A quantizer  assigns  each  sample  to  a particular 
interval  according  to  its  amplitude  (quantization  is 
sometimes  considered  to  be  sampling  in  amplitude).  For  an 
individual  sample,  the  interval  is  a portion  of  the  real 
line  (see  fig.  3-la).  However  if  two  samples  at  a time  are 
considered,  then  a quantizer  assigns  pairs  of  samples  to 
regions  in  two-space,  as  shown  in  fig.  3-lb.  Similarly,  a 


quantizer 

assigns 

vectors  of 

N samples 

to  reqions  in 

N-space , 

RN.  If 

the  vector 

components 

are  quantized 

independently,  the  resultant  region  in  RN  is  rectangular. 
An  estimate  of  the  quantizer  input  based  on  its 
corresponding  output  region  would  be  one  of  the  points 
within  that  region.  The  goal  in  this  chapter  is  to  find 
the  estimation  point  which  is  most  similar  to  the  quantizer 
input,  with  respect  to  a given  error  criterion.  It  should 
be  mentioned  that  no  attempt  is  made  to  find  optimum 
quantization  regions.  Rather,  a method  is  derived  for 
finding  the  optimum  restoration  point  within  an  arbitrary, 
given  region.  The  restoration  concepts  can  best  be 
understood  in  the  context  of  the  following  simple  example. 


3.2  Vector  Restoration  Example 

Consider  a two  level  (one  bit)  quantization  of  a 
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sampled  signal  limited  in  amplitude  to  a finite  range. 
Specif icallYf  for  a signal  limited  to  the  range  [-1,1],  let 
the  output  of  a guantizer  be  the  information  that  the 
signal  is  either  in  [0,1]  or  in  [-1,0).  Next  assume  that 
two  successive  outputs  at  times  t^  and  t 2 are  as  shown  in 
fig.  3-2a. 

x ( tx ) e [0,i] 

x ( 1 2 ) e 1-1,0) 

If  nothing  else  is  known  about  these  signals,  then  the  best 
restoration,  as  fig.  3-2b  shows,  would  be 

y(t1)=i/2 

y(t2)=-l/2 

However,  if  the  signals  are  samples  from  a random  process 

with  a known  probability  distribution,  p(x) , such  as  the 

one  in  fig.  3-2c,  then  a better  restoration  would  be  the 

mean  values  of  each  interval  according  to  this  distribution 

(see  fig.  3-2d).  (The  restoration  in  fig.  3-2b  would 

correspond  to  a uniform  distribution.)  Finally,  if  the 

signals  at  t and  t are  known  o be  correlated,  then  the 
12 

restoration  can  be  further  improved  by  utilizing  this 
correlation.  Figures  3-2e  and  3-2f  show  the  restoration 
points  for  positive  correlation  and  negative  correlation, 
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respectively. 


3.3  Nonlinear  Estimator 

To  characterize  the  above  concepts  mathematically,  an 
expression  for  the  quantization  error  must  first  be 
derived.  Let  a vector  of  N samples  of  a random  process  be 
denoted  as  x,  where 


x=(xx  ,x2  ,. . . ,xN)  (3.1) 

and  assume  that  this  vector  is  statistically  described  by 

its  probability  density  furction,  p(x) . If  this  vector  is 

quantized,  according  to  the  techniques  of  Sec.  3.1,  then  it 

is  assigned  to  one  of  M regions  in  RN , denoted  as  , for 

m=l,2,...,M.  Next  assume  that  the  quantizer's  desired 

output  is  its  input  signal.  This  is  a reasonable 

assumption,  for  in  most  communication  systems  the  quantized 

signal  is  intended  to  be  an  instantaneous  replica  of  the 

input  signal.  Thus,  it  can  also  be  assumed  that  x € R , 

— m 

i.e.,  that  the  input  vector  x is  assigned  to  the  region  in 
which  it  is  contained.  Note  that  these  assumptions  are 
made  for  convenience  only  and  place  no  restrictions 
whatsoever  on  the  subsequent  solution. 

The  point  within  R that  is  chosen  as  the  restored 

m 

value  of  the  quantizer  input  is  denoted  as  where  ^ is 
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a vector  since  it  represents  a point  in  N-space.  The 


instantaneous  error  caused  by  restoring  a quantized  signal 


to  the  point  ^ is  then  x-^.  The  restoration  performance 
is  to  be  determined  by  evaluating  a function  of  this  error, 


e(*)*  The  goal  is  to  minimize  the  mean  value  of  this  error 


function,  for  all  possible  quantizer  inputs,  by  the  proper 


choice  of  restoration  points,  The  expression  for  the 


total  error  which  must  be  minimized  is  then 


M f 

■-£  /e(x- 

m=l  J 


p(x)dx 


(3.2) 


The  error  weighting  function,  e,  is  arbitrary,  but  it  is 


usually  a nonnegative  function  because  instantaneous  errors 


of  opposite  sign  should  not  cancel  each  other. 


Many  researchers  have  proposed  solutions  which 


minimize  eq.  3.2  for  different  choices  of  the  error 


weighting  function,  but  all  of  their  results  are  predicated 


on  x being  a one-dimensional  or  scalar  input.  Of  these 


results,  the  most  significant  ones  have  been  obtained  by 


Bruce  [1]  who  used  dynamic  programing  techniques  to  find 


the  optimum  regions,  , and  restoration  points,  for 


arbitrary  error  weighting  functions.  His  results 


generalized  those  obtained  earlier  by  Max  [2]  for  a squared 


error  weighting. 


yi 


'■*  -A  * 


utilized  here. 


The 


A mean-square  error  criterion  is  also 
functional  form  of  this  error  weighting 


is 


T- 


e(2L-Ym)=Tr{{-"^)  (-"^)  1 


(3.3) 


where  Tr { • } denotes  the  trace  of  a matrix 
function  e(x-i^)  is  now  a monotonic  function 
a simple  assumption  on  p(x)  (that  p(x)  is 


The  error 
Then,  under 
not  entirely 


discrete) , the  error  expression 
by  the  techniques  of  differentia 
assumption  about  p(x),  the  re 
surface  $ can  be  found  from 


in  eq.  3.2  can  be  minimized 
1 calculus.  By  making  this 
lative  extrema  of  the  error 


(3-4) 


for  m=l , 2 , . . • i M . 
eq.  3.4  yields 


Substituting  equations  3.2  and  3.3  into 


(3.5) 


for  m=l,2, 


,M.  This  equation  can  be  rearranged  and  the 
, solved  for  to  obtain 


restoration  point,  y 


xp(x)dx 


V 


m 


(3.6) 


m 

for  m=l,2,...,M.  Equation  3.6  is  thus  a minimum 
mean-square  error  estimate  of  the  input  to  a quantizer, 
based  on  the  quantizer  output  and  the  a priori  statistics 
of  the  input.  Closer  inspection  of  this  equation  reveals 
that 


y =E { x | x € R } 
Tn  - - m 


(3.7) 


for  m=l , 2 , . . . ,M.  That  is,  ^ is  the  conditional  mean 

estimate  of  x conditioned  on  the  nonlinear  information  that 

x has  been  quantized  to  the  region  R . Because  of  this 
— m 

conditioning,  represents  a nonlinear  res'  oration  of  the 

quantized  signal.  The  restoration  is  optimum  only  with 
respect  to  a mean-square  error  criterion,  but  it  is 
completely  independent  of  the  specific  form  of  the 
quantizer . 

Note  that  eq.  3.6  requires  that  a multidimensional 
density  function,  p(x),  be  provided  in  order  to  solve  for 
the  optimal  restoration.  Since  a gaussian  density  is  the 
only  known  multidimensional  density  for  correlated  signals, 
only  quantized  gaussian  random  processes  might  be  optimally 
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restored.  If  other  multidimensional  densities  can  be 
derived,  then  systems  which  are  based  on  these  densities 
miqht  also  be  optimized.  The  subsequent  chapter  considers 
this  derivation  possibility  and  presents  a technique  for 
generating  correlated  multidimensional  densities.  The 
results  are  then  utilized  in  succeeding  chapters  to  analyze 
several  common  communication  systems. 

3.4  Estimation  Covariance 

The  performance  of  the  estimator  derived  in  the 
previous  section  can  be  determined  by  computing  the 
estimation  error  covariance  matrix.  The  covariance  of  the 
estimator,  based  on  the  quantized  measurement  information, 
can  be  found  from 


cov(x  l X e V=E{(*-V  (2!-VTl  * e v (3.8) 


By  comparison  with  equations  3.6  and  3.7  thib  can  be 
expressed  as 


cov(x  | x 6 R ) 
— — m 


T 


p (x ) dx 


p (x ) dx 


(3.9) 


The  numerator  can  then 
equation  simplified  to 


be  expanded  and 
obtain  the  final 


the  resulting 
form  for  the 
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estimation  covariance 


cov(x  | x 6 R )=- 
— — m 


/, 


xx  p(x)dx 


R 


m 


/ p(x)dx 

Jr 


T 

*m*m 


(3.10) 


m 


Using  this  result,  the  average  quanti7.s.cion  error  remaining 
after  the  minimum  mean-square  error  restoration  of  Sec.  3.3 
is 


M 


£ = Y'  Tr  { cov  (x  | x 6 R )}Pr(x  0 R ) 
m=l  - - m - m 


(3.11) 


where  Pr (A)  denotes  the  probability  of  occurrence  of  the 
event  A ana 


Pr (x  0 R ) 
~ m 


= / p ( x ) dx 

Jr 


(3.12) 


m 


Then  using  eq.  3.10,  the  error  becomes 


M 


£=Tr  £ 

*m=l 


/, 


xxTp(/)dx  - y_y^Pr(x  0 R ) 
~ ■*TndTn  — m 


m 


(3.13) 


which  simplifies  to 


M 


<?=Tr{C  0 RJ} 


m=l 


m 


(3.14) 
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where  the  matrix  C is  the  covariance  of  the  input  vector 
x.  Alternatively,  this  same  result  could  have  been  derived 
by  substituting  equations  3.6  and  3.3  into  eq.  3.2  and  then 
simplifying  the  result.  The  error  expression  in  eq.  3.13 
is  used  later  to  evaluate  the  performances  of  various 
optimal  restoration  techniques  that  are  based  on  quantized 
measurements . 
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CHAPTER  4 


MULTIVARIATE  PROBABILITY  DENSITY  FUNCTIONS 


4.1  Introduction 


Multivariate  distributions  of  dependent  random 
variables  are  often  required  in  communications  problems  for 
the  analysis  of  stationary  random  processes.  The 
distributions  arise  whenever  a continuous  random  process  is 
sampled — a fundamental  characteristic  of  all  digital 
systems.  The  probability  density  function  of  N samples  of 
a random  process  X(t)  is  usually  denoted  as 


p(X  ,X,..., X ) 

1 2.  N 

This  expression  is  very  general  and  represents  an  infinite 
number  of  possible  density  functions,  only  a few  of  which 
have  ever  been  found.  In  fact,  for  N>2  and  correlated 
variables,  there  is  only  the  familiar  multivariate  gaussian 
distribution.  Fortunately  the  gaussian  is  the  most  useful 
distribution  (by  virtue  of  the  Central  Limit  Theorem  making 
it  the  limiting  distribution  for  many  additive  processes). 
There  are,  however,  many  processes  which  it  does  not  model 
well.  For  these,  the  lack  of  other  known  distributions 
often  means  that  a multivariate  gaussian  is  used  by 
default,  or  that  the  dependence  between  the  variables  is 
ignored  and  the  density  is  written  as  the  product  of 


.7-i  7 


. Iiupp 


independent 

marg inal 

densities.  In 

either 

case 

inaccuracies 

result . 

Thus,  this  chapter 

addresses 

the 

problem  of 

finding 

the  multivariate 

densities 

of 

non- gaussian 

variables 

whose  marginal 

densi t ies 

and 

correlation  function  are  known. 

4.2  Characteristics  of  Multivariate  Distributions 

The  probability  density  function  of  a single  random 
variable  is  uniquely  determined  by  specifying  all  M+l  of 
its  moments,  l£M<^o  , denoted  by  [1] 

E { xm } 

for  m=0 , 1 , 2 , . . . ,M.  Uniqueness  holds  under  the  following 
easily  satisfied  conditions:  when  all  of  the  moments  are 

finite  and  when  a power  series  with  the  moments  as  its 
coefficients  converges.  To  uniquely  specify  the 
multivariate  density  of  N random  variables  requires  (M+1)N 
joint  moments,  similarly  denoted  as 

E{x™l  x™2  . ..xJJn  } 

for  m^  ,m2  , . . . ,mN=0 ,1 , . . . ,M.  For  a g..ven  random  process  it 
is  usually  possible  to  determine  the  distribution,  p(x^) , 
of  a single  random  variable  (a  marginal  density  function  of 
the  multivariate  density) , and  to  measure  its  correlation 
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with  the  other  similarly  distributed  variables.  If  there 

are  N variables,  then  the  N marginal  densities  provide 

2 

N ( M+ 1 ) moments  and  the  correlations  provide  N -N  more.  The 

remaining  mixea  moments  are  undetermined.  If  these  are 

chosen  arbitrarily,  then  there  are  an  infinite  number  of 

2 

possible  multivariate  densities  which  have  the  same  NM+N 
moments  specified  above.  However,  these  densities  differ 
from  each  other  only  in  their  higher  moments,  and  thus  are 
very  similar. 

A valid  probability  density  function  must  satisfy  the 
following  six  conditions  [2]: 


Property  1. 


P(x^  , * • • , (4.1) 


Property  2.  /.../  p(x.  ,x_  , . . . ,x„)dx.  dx_  . . .dx  =1  (4.2) 

* I 2 N 1 2 N 


Property  3.  / p ( x , x , . . . , x ) dx . =p (x  , . . . , x . fx.,f...,x„) 

/ 1 2 N 1 1 1-1  1+1  N 


(4.3) 


for  i=l , 2 , . . . ,N. 


Property  4.  To  accurately  represent  a given  random  process 
with  correlation  matrix  C,  the  density  function  must  be 
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able  to  generate  the  elements  of  C from  the  equation 


^-I-h 


• p(x^  ,*2  i • • • r )dx^dx2»«»dXj(j  (4.4) 


Property  5.  When  a stationary  random  process  X(t)  is 
sampled,  the  samples,  in  general,  are  correlated  with  a 
correlation  matrix  C.  If  the  sample  spacing  becomes  so 
large,  however,  that  the  samples  are  uncorrelated  and  C is 
a diagonal  matrix  (assuming  X(t)  has  no  strictly  periodical 
components) , then  the  samples  are  also  independent.  Then 
the  density  function  should  be  factorable  into  a product  of 
independent  density  functions,  as 


P ( x^  , > •••  i xj^  ) — P ( x^  ) p ( *2  ) ♦ • • P ( ) 


(4.5) 


Property  6.  If  two  of  the  random  variables  are  identical, 

e.g.  x.=x.,  then  the  density  must  become 
i 3 


p(xx  ,x2  ,.. .,xN)*8(xi-x  )p(xx  ,. . . ,xi_jxi+f • • • fxN)  (4.6) 


4.3  Prior  Multivariate  Density  Research 


A few  specialized  multivariate  densities  have  been 
found  in  the  past.  In  1945  S.  0.  Rice  [3]  derived 


two-dimensional 


Rayleigh  and  sine-cosine  densities. 
Subseauent  researchers,,  such  as  Gumbel  [4]  and  Parzen  [5], 
found  multivariate  densities  based  on  known  marginals,  but 
their  densities  failed  to  satisfy  Properties  4 and  6 and  so 
could  not  represent  a given  random  process.  The  most 
significant  results  to  date  were  obtained  by  Beckmann  [6] 
who  developed  a general  technique  for  constructing  a 
two-dimensional  density  from  its  marqinals  and  correlation 
coefficient.  Beckmann's  method  required  each  marginal 
density  to  be  the  weighting  function  of  a system  of 
orthogonal  polynomials.  Unfortunately  this  resulted  in  a 
joint  density  containing  an  infinite  series  that  often  did 
not  have  a closed  form  solution.  Also,  by  the  inherent 
nature  of  orthogonal  polynomials,  the  results  were  not 
extendable  to  more  than  two  dimensions. 

4.4  Derivation  of  Multidimensional  Densities 

To  derive  densities  in  higher  dimensions,  consider 
first  an  N-d imens ional  gaussian  density.  Assume,  without 
loss  of  generality,  that  its  mean  is  zero  and  its 
covariance  is  C.  Then  this  density  can  be  written 

p(x)  = ( 2tt ) _N/  ref  1/2exp(-|xTC-1x)  (4.7) 

It  is  possible,  by  a suitable  linear  transformation,  to 
express  this  density  in  a new  coordinate  system  in  which 
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the  components  are  statistically  independent.  The  density 


transformation  is  effected  through 


pv  (y)=pY(£> 


(4.8) 


T 

y=Ex  x 


(4.9) 


where  E is  the  solution  to  the  eigenvalue  equation 


EA=CE 


(4.10) 


with  A a diagonal  matrix  of  eigenvalues.  Since  C ij 
symmetric  and  a covariance  matrix,  the  eigenvalues  are  real 


and  nonnegative.  Also  the  eigenvectors  which  comprise  the 


columns  of  E can  be  chosen  to  be  orthonormal,  i.e. 


ete=i 


(4.11) 


Now  substituting  eq.  4.9  into  equations  4.8  and  4.7  yields 


p (y)  =p  (Ey)  |E_1  | = (2tt)"N/2  I E|  1 | C |'1/2  exp  ( -^yT  ETC  1 Ey ) 
Y X 


■ft  ’>  v>'  ■ . 

..*.'** a 


(V)  = (27t)_N/2  | Af1/2exp(-jXT  A_1  1)  (4.12) 

Notice  that  i£  A in  this  last  equation  is  an  identity 
matrix,  then  the  original  covariance  matrix  C must  also  be 
an  identity  matrix,  for 

C = EAET  = EET  = I (4.13) 


The  diagonal  elements  of  A are  just  the  variances  of 
the  components  in  the  new  coordinate  system.  If  the 
components  are  scaled  by  the  square  root  of  these 
variances,  i.e.  = y^x/ X ^ , then  each  new  component  has 

unit  variance.  Hence  rotating  these  scaled  coordinates  by 
the  inverse  of  the  original  transformation 


yields  an  uncorrelated  unit  variance  probability  density 
function 


1(w)  = (27r)"N/2exp(4wTEA1/2A"1A1/2E"1w) 


= (2  7r)_N//2exp(-^wTw) 


(4.15) 


Reversing  the  above  procedure  thus  leads  to  a method  for 

generating  a correlated  multidimensional  density  function 
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from  an  uncorrelated  one. 


Reiterating  the  steps  above  in  proper  order  yields  the 
following  general  technique: 

(1)  Write  an  uncorrelated  N-dimensional  density  as  a 
product  of  its  known  marginals 

p(x)=p(x1 ) p(x2 ) . . .p(xN ) (4.16) 

(2)  Define  E such  that  E is  orthonormal  and 

ETCE=A  (4.17) 

where  A is  diagonal  and  C is  the  desired  covariance  matrix. 
The  matrices  E and  A are  respectively  the  matrix  of 
eigenvectors  and  the  matrix  of  eigenvalues  of  C. 

(3)  Transform  x by  E , i.e.,  y=E  x=E  x,  so  that 

PY(¥)=lE|px(Ex)=px<Ex)  (4.18) 

(4)  Scale  the  components  of  ^ by  the  square  roots  of  the 
eigenvalues;  then 


z=A 


(4.19) 
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and 

-1 

Pz  (z)  = |A!"1/2EV  (^1/2  z)  = lAf1/2px  (EA“1/2  z) 

(4.20) 

(5) 

Inverse  transform  z by  w=Ez  to  obtain 

% (w)  = ] E |Tpz  (ETw)  = |A|'1/2px  (EA'1/2ETw) 

% (wMAlp^  (Aw) 


where 


a=ea'1//2  et 


(4.21) 


(4.22) 


To  illustrate  the  effect  of  each  step  of  the  procedure, 
fig.  4-1  shows  contour  plots  of  two-dimensional  gaussian 
and  laplacian  densities. 

4.5  Proof  of  Density  Properties 

The  function  pw (w)  in  eg.  4.21  has  correlation  matrix 
C and  satisfies  the  conditions  for  a valid  probability 
density  function  as  is  next  shown. 


GAUSSIAN  DENSITY  LAPLACIAN  DENSITY 


Uncorrelated  Density  Contours 


Coordinate  Rotation 


Eigenvalue  Scaling 


Inverse  Coordinate  Rotation 


Figure  4-1.  Contour  plots  of  the  two-dimensional  density 
transformations  used  to  generate  correlated  densities. 


Property  1. 
for  all  x e 

Prooerty  2. 


Property  3. 
Property  4. 


Pioperty  5 

Property  6 

lim  a . . = 
ii 


pw(w)  is  clearly  nonnegative  because  Px(x)> 


RN  and  Aw  € RN. 


/pw(w)dw=  / I A | p^ 

,n  Jja 


w 


(Aw) dw 


/Px  (x)d* 


= 1 


J pw  dVpNi  '*•*  'Viwi+i  • • • 'V 


E{wiw.  }=E{  (A~1x)i  (A  1x)j) 


- 1 , T,  -1 
=A^  E { xx  } Aj 


( EA1/2  ET  EA1/2  ET ) i;j 


-(£»u 


C=I  implies  that  A=I^;  then 


(4.23 


(4.24 


( 4 . 21 


I A|px  (Aw)  =PX  (w)=p(w1  )p(w2)  . . .p(wN)  (4.2 


. IE  x.=xj.  then  j-aj . . H 


a , . = oo  . Hence  lim  |A|  = oo  and 
ID 


v 


V 


Thus 


0 ; X.  >X  . 
1 3 


lim  (w)  = J oo  ; xi=x;. 

0 ; X <x 


Xi~*  Xj 


= 5 { x^  — x . ) 


(4.28) 


Unfortunately,  Property  3 is  not  satisfied,  i.e.,  the 
correlated  density  function  no  longer  has  the  desired 
marginal  distributions.  In  fact,  its  marginal 
distributions  are  a function  of  the  correlation  matrix  C. 
However,  as  is  shown  in  more  detail  later,  the  correlated 
density  function  found  above  remains  a good  representation 
for  a given  random  process,  and  is  a good  approximation  to 
the  exact  (but  unknown)  probability  density  function. 

4.6  Examples  of  Multidimensional  Densities 

To  demonstrate  the  utility  of  the  derivation 
procedure,  some  examples  are  now  presented. 


* jm 


i 


A.  Gaussian  probability  density 


Px  ( x ) = (2  tt)-N//2  exp(-|xTx) 


(4.29) 


(w)  = I A l_  1/2  (27r)“N/2exp(-jWTATAw) 


= |C|"1/2  (2T)“N/2exp(-jWTEA_1/2  ETEA_1/2ETw) 


ICf 1/2  (2*)"  N/2exp(-^wTC-1w) 


B.  Laplacian  probability  density 


N 


PX(i)*^72eXP(-VJi?1IXil) 


p„(!!i-|£|'1/2-572exp(-v5.Ei 


N 


a. , w 


k-1  lk  k 


C.  Rayleigh  probability  density 


px<*>= 


x)  = j^J]^xiU(xi  )j  exp(-^-xTx) 


(4.30) 


(4.31) 


(4.32) 


(4.33) 


rW 


( w) = i c r1/2  ( x>ik  wk ) u < k5xaik  \ >]  exp ( ~hT  £ ^ 


w)  (4.34) 


D.  Maxwell  probability  density 


Px(x)  = (2/ir) 


N/2 


N 


U(x.) 


1 T 
exp ( x x) 


(4.35) 


( 2/t  ) 


N/21 


N 


N 


TT ( £ aikwk ) 

1=1  k=l  k k 


, N 

U(  Z aikwk) 

k=l  lk  K 


exp  ( -WC\)  (4.36) 


Figures  4-2 
gaussian , 


through  4-9  contain  plots  of  two-dimensional 
laplacian , and  Rayleigh  densities  for  both 


uncorrelated  and  correlated  variables. 


4.7  Marginal  Densities  and  Random  Process  Simulations 


As  stated  previously  (see  Section  4.5),  the  correlated 
multidimensional  density  function  in  eo . 4.21  and  the 
uncorrelated  density  in  eg.  4.16  have  differtnt  marginal 
density  functions.  In  generating  the  correlated  density  by 
steps  (1)  through  (5),  the  marginal  density  becomes  a 
function  of  the  correlation  matrix  C.  For  example,  a 
two-dimensional  correlated  laplacian  density  can  be  written 


1 

P(x  ,y)  = f=yexP 

2s  s \ 1-r^ 
x y 


1 

yj2(l-r2  ) 


hy. 
s y 


ay. 
s „ 


(4.37) 


where 
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p(X,  Y) 


FIGURE  4-4. 


TWO-DIMENSIONAL  GAUSSIAN  DENSITY 
Correlation  - -0.  8 


f .) 


% 


sx=E{y2} 


(4.38) 


Sy=E{y2} 

E{  xy} 

r = 


a=Vl+r  + y/l-r 


b=Vl+r  - vT-r 


Its  marginal  distribution  is 


(4.39) 


(4.40) 


(4.41) 


(4.42) 


p(x)  = 


2s  \/2  ( 1-r  2 ) 


( a | exp(-2^|£-  I j _ | b|  exp  (-2^ 


_2L 

b 


(4.43) 


which  is  a function  of  the  correlation,  r.  It  becomes  the 
classical  laplacian  distribution  only  when  the  correlation 
is  zero,  i.e. 


P(x) 


r=0 


V2  \ x | 


(4.44) 


However,  non-zero  correlations  in  ea.  4.43  result  in 
marginal  distributions  that  are  very  similar  to  the  one 
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described  by  eq.  4.44.  Fig.  4-10  shows  these  marginal 
distributions  for  various  correlation  coe f f icients— the 
differences  are  minimal  for  even  a larqe  change  in  the 
correlation.  Thus  either  function  could  adequately  model  a 
random  process  such  as  the  one  characterized  by  fig.  4-11. 
This  figure  represents  the  measured  distribution  at  the 
output  of  a DPCM  coder  having  a sampled  image  as  an  input. 
This  particular  image  has  a measured  average  correlation  of 
0.4  between  successive  DPCM  samples,  resulting  in  the 
experimental  two-dimensional  distribution  shown  in 
fig.  4-12.  Using  r=0.4  in  eq.  4.37,  this  distribution  can 
be  modeled  as  a two-dimensional  laplacian  density.  This 
density  is  plotted  in  fig.  4-13. 


Similarly,  fig.  4-14  shows  the  marginal  distributions 
of  a two-dimensional  Rayleigh  density  for  several  different 
correlations.  The  curves  are  plotted  from 


P ( x ) = * exp 
2s 

x 


r 2, a 2~ 

2 

X 

“eXP 

x (6— r) 

2 

2sx 

! 

X 

2sx(l-r  )J 

(l-Vl-r2) 


+§-exp 

x 


2 .b  ,2 

* (r-r) 

? 2 
2sJ(l-r  ) 


(l+yji-t7) 


r-  2 
hr  .x 

+rVr<-- 


1) 


erf 


, a » 
*<E“r) 

-er  f 

,b 

x(r- r) 

— 

sx\/2(l-r2)_ 

U(x)  (4.45) 
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Figure  4-10.  Marginal  distributions  of  a correlated 
two-dimensional  laplacian  density  function  for  various 
correlation  coefficients  (r)  . 


FIGURE  4-13.  TWO-DIMENSIONAL  LAPLACLAN  DENSITY 
FUNCTION  USED  TO  MODEL  THE  TWO-DIMENSIONAL  DPCM 
SIGNAL  SHOWN  IN  FIG.  4-1Z;  CORRELATION  =0.4 
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T--N  - > 


which  is  found  by  integratina  the  two-dimensional  Ravleiqh 
density 


P(x ,y) = 


,ax  by  t <ay  bxwl/ax_bY,nna^_bX' 

4s  s (1-rV75  s x"s  y s Y S x SxSy  SyS- 
x y 


•exp 


2 _ 2 
1 2r  xy.  + JL-) 

2^ ( 7 2 

2 t 1-r  ) s s s s 
' x x y y 


(4.46) 


The  parameters  s , s , r,  a,  and 

x y 

previously  defined  in  equations  4.38 
when  the  correlation  between  x and 
marginal  density  reduces  to 


b are  the  same  as 
to  4.42.  As  expected, 
y becomes  zero,  the 


P(x) 


2 

-S-tJIHx) 

2s  * 
x 


(4.47) 


the  familiar  one-dimensional  Rayleigh  density.  However  the 
figure  shows  that  the  shape  of  the  marginal  density 
function  is  insensitive  to  changes  in  the  correlation. 
Hence  a random  variable  approximately  described  by  eq.  4.47 
could  be  characterized  equally  well  (though  not  as 
concisely)  by  eq.  4.45.  Similar  results  are  obtained  for 
marginal  distributions  of  the  higher  dimensional  densities 
described  by  eq . 4.34.  This  equation  therefore  provides  an 
effective  model  for  a multidimensional  correlated  Rayleigh 
random  process. 
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Thus,  the  technique  developed  above  generates,  from  a 
given  one-dimensional  density,  a multidimensional  density 
that  possesses  a desired  correlation  function.  Because  the 
one-dimensional  density  function  and  correlation  properties 
of  a random  process  are  usually  measurable  and  known,  the 
multidimensional  density  can  be  derived  from  this 
information  and  utilized  to  model  the  random  process.  This 
model  can  then  be  employed  in  the  analysis  and  simulation 
of  digital  systems  which  operate  on  the  random  process. 
The  potential  applications  for  this  modeling  are  widespread 
and  several  are  analyzed  in  subsequent  chapters. 


REFERENCES 


1.  A.  Papoulis,  Probability,  Random  Variables,  and 
Stochastic  Processes,  McGraw-Hill,  Inc.,  New  York,  1965, 
pp.  157-158. 


2.  Ibid.,  chap.  8. 


3.  S.  0 Rice,  "Mathematical  Analysis  of  Random  Noise," 
Bell  System  Technical  Journal,  vol . 24,  1945. 


4.  E.  J.  Gumbel , 


"Distributions  a 
dont  les  marges  sont  donnees," 

1 ' Academie  des  Sciences,  Paris, 

pp.  2717-2720 . 


plusieurs  variables, 
Comptes  Rendues  de 
vol.  246,  1958, 


5.  E.  Parzen,  Modern  Probability  Theory  and Its 

Applications,  Wiley,  New  York,  1960,  p.  279. 

6.  P.  Beckmann,  Orthogonal  Polynomials  for  Engineers — and 
Physicists,  Golem  Press,  Colorado,  1973,  pp.  163-187. 


69 


CHAPTER  5 


QUANTIZATION  AND  RESTORATION  OF  GAU  ^ I AN  SAMPLES 

In  Chapter  3 , an  estimation  equation  was  derived  for 
the  minimum  mean-square  error  restoration  of  quantized 
samples  from  a random  process.  Solution  of  the  estimation 
equation  is  dependent  upon  knowledge  of  both  the  structure 
of  the  quantizer  and  the  multidimensional  probability 
density  function  of  the  sampled  random  process.  For  many 
communication  systems,  the  underlying  random  process  can  be 
described  by  a gaussian  probability  density  function.  A 
gaussian  density  arises  naturally  in  many  applications 
because  it  is  the  limiting  distribution  for  all  additive 
random  processes,  and  the  addition  of  random  variables 
occurs  often  in  communication  systems.  A gaussian  random 
process  is  thus  assumed  to  be  the  input  to  a quantizer  and 
an  optimal  restoration  is  then  derived  based  on  this 
assumption . 

5.1  Estimation  of  Quantized  Gaussian  Samples 

Repeating  one  of  the  major  results  of  Chapter  3,  the 

restoration  of  a vector  of  samples  x,  which  has  been 

quantized  to  one  of  M regions,  R , can  be  found  from 

m 

k 


i 


xp ( x ) dx 


I o(x)dx 

V " 


(3.6) 


for  m=l,2,...,M,  where  p(x)  is  the  probability  density 


function  of  the  sampled  random  process.  For  a gaussian 


random  process,  p(x)  can  be  written 


P(x)  = (2t)“N/2  | Cf 1//2  exp(-^xT C 1x) 


(5.1) 


where  C is  the  covariance  matrix  of  x and  the  mean  is 


assumed  tc  be  zero.  If  the  region  R is  rectangular,  i.e., 


if  each  component  of  x is  quantized  individually,  then  the 


region  can  be  expressed  as 


R ={x  | x,6  [a  . ,b  .],  i-1 , 2 , . . . ,N} 
m - i mi  mi 


(5.2) 


for  m=l , 2 , . . . ,M.  For  notational  simplicity  only  a single 


region  is  henceforth  considered  and  the  subscript  m is 


dropped,  leaving 


R= { x | x € [a,b|  ) 


(5.3) 


Substituting  eq.  5.1  into  ea.  3.6  yields 


T 

' ,v> 


(5.4) 


"x(2t r)  N//2  | C r1//2exp(-ixTC_1x)dx 
a 

~ ( 2 7T)  “N/ 2 I cr1//2exp(-^xTC_1x)dx 
a 


Unfortunately,  no  known  analytical  solution  exists  for 
either  of  the  integrals  in  this  eauat ion . Curry  [1)  has 
obtained  an  approximate  solution  for  finely  quantized 
values  of  x^,  i.e.,  for 


b -a . < s . (5.5) 

ii  l 


for  i = l , 2 , . . . , N , where  si  is  the  standard  deviation  of  the 
( i) th  component  of  x.  His  approach  is  to  approximate  the 
gaussian  density  by  the  first  three  terms  of  its  Taylor 
series  expansion  about  the  midpoint  of  the  region  R.  The 
integrations  can  then  be  performed,  with  the  result  that 


E{x  I x 6 R}=( I-DC_1) (b+a)/2 


(5.6) 


where  the  diagonal  matrix 


D= 


(bi-ai} 


12 


5.  . 

ID 


for  i,]-l,2,...,N. 


(5.7) 
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An  exact  solution  to  eq.  5.4  can  be  obtained  when  the 
components  of  x are  uncorrelated.  In  this  case  the 
covariance  matrix  of  x can  be  expressed  as 


C=  { s.  8.  . } 


(5.8) 


for  i , j = l , 2 , . . . ,N , so  that  the  restoration  equation  can  be 
rewritten  as 


y=- 


f xexp(-i  ^ x?/s?)dx 
; a-  \ 2 i = 1 1 V - 


II  exp(i?i*Hd- 


(5.9) 


Each  multidimensional  integral  can  then  be  separated  into  a 
product  of  one-dimensional  integrals,  so  that  common 
factors  can  be  cancelled,  leaving 


1 


Considering  the  ( i ) t h component  of  this  vector,  since  all 
of  the  components  are  identical  in  form,  the  integrals  can 
be  evaluated,  as  <*> 


<*>  In  accordance  with  the  usual  convention, 


ft 


A similar  computation  without  the  zero-mean  assumption 
shows  that,  if  E{x}=u,  then 


\yZ2sJ  \V2sJ 


(5.12) 


N 

erf/^^Verffll^A 

\V2sJ  Us  J 


(b  -u  .) 


Gaussian  variables  which  have  been  decorrelated  by  means  of 
a Karhunen-Loeve  transformation  [2]  and  then  quantized  can 
be  restored,  according  to  a minimum  mean-square  error 
criterion,  by  utilizing  these  last  two  equations. 


also  exists 


An  exact  analytical  solution  to  eq.  5.4 
when  an  estimate  of  a single  vector  component,  x^ , is 
desired  based  upon  two  types  of  information: 

1.  the  other  components,  x2  ,x2  ,. . . ,xN_j,  which  are 
known  completely  (quantized  with  an  infinite 
number  of  bits) 

2.  the  quantizer  output,  which  nonlinearly  specifies 

the  interval  containing  x., 

N 

To  derive  this  solution,  consider 


t exp(-±t  C t) dx 


exp(-|tTC_1 1)  dx 


t exp(-^tTC-1t)  dx 
2 — — — N 


exp  ( “tTC- 1 1)  dx 
2—  — — N 


Nx  exp(-'S‘tTC-1t)  dx 
N 2 N 


Nexp(-ltTC  -*-t)dx 


2 


Denoting  the  elements  of  C_  by 


11  12  ‘ ‘ 


the  exponential  term  the  Nth  component  of  eg 

becomes 


=exp 


N-l 


N-l 


“l  ;SlirlD  aD  +a2  ^_lr2  j “j  + 


N-l 

+aN-l  £Vl,j 


•exp 


”2  *N  (rjN+rNj  )+rNN*N  ] 


(5.17) 


The  first  exponential  in  eq.  5.17  is  a constant  and  is 
common  to  both  the  numerator  and  denominator  of  the  Nth 
component  of  y in  eq.  5.15.  Hence  this  factor  can  be 
cancelled,  leaving 


J 

^bN 

vxp 

1 

~2 

(rjN+rNj  )+rNN*N 

dXN 

yN  = " 

v> 

r b 

f NexP 

' % 

-s-[ 

N-i  21 

*Nj^iaj  ^rjN+rNj  ^+CNNXN  J 

dXN 

(5.18) 


Completing  the  square  within  the  exponential, 
exponential  becomes 


each 
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1 


r 


exp 


i 

I CNN  ("  2 XN  ^ / ! N - 1 \ 

MVr^  :?i‘V>+%,+(2— N £,VW) 


/ 1 N-l  v 2 

"fe  -j^VW 


=exp 


rMM  r i n-i  "i 

2 LN  2rNN  i?i^  r^N  'N^  . 


>exp 


rMM  r n N-l 

■Sfe  &v>*v] 


+ r . ) 
]N  N] 


(5.19) 


Aqain,  the  last  exponential  is  constant  and  is  common  to 
both  the  numerator  and  denominator.  Hence  it  can  be 
cancelled.  Next,  performing  the  substitutions 


N-l 


* Vv 


N 2r«N  fr  i 3N 


(5.20a) 


du=dx 


N 


(5.20b) 


T" 


(5.23) 


P ~~2^(bN  + rf)  _exP  __7i(aN+w) 

erfV^"(bN+W)  "erf  (aN  <-w) 


where 


N-l 

w=-^—  T a.r  . 
rNN  ] = 1 3 N3 


(5.24) 


If  XN  is  quantized  to  an  infinite  number  of  bits 


(a-b),  then 
N N 


Y =a 
N N 


(5.25) 


as  is  expected.  If  x is  quantized  to  zero  bits,  it: 


quantization  interval  is  the  real  line  (-^=^  = 00 ) r and 


then  its  estimate,  y , is 

N 


. N-l 

y =— — £ a.  r . 
r NN  j=l  3 N3 


(5.26) 


‘ 


This  result  is  identical  to  one  obtained  by  Pratt  [3]  in 
estimating  an  unknown  spectral  value  based  on  known 
spectral  components.  (The  equivalence  between  these 
results  is  shown  in  the  next  section.)  However,  ea.  5.23 
is  a more  general  result  in  that  it  can  be  utilized  to 
estimate  components  that  have  been  quantized  to  any  number 
of  bits  by  an  arbitrary  quantization  scheme. 

By  comparing  the  Nth  component,  y^,  of  the  zero-mean 

version  of  eq.  5.12  (the  solution  for  uncorrelated 

variables)  with  y found  from  eq.  5.23  (the  solution  for 
N 

correlated  variables),  it  can  be  seen  that  the  correlation 

effects  are  embodied  in  the  variable,  w.  In  fact,  y^  from 

eq.  5.23  reduces  to  y„  from  eq.  5.12  when  w=0,  which  in 

turn  occurs  when  the  correlation  is  zero.  The  variable  w 

is  simply  the  negative  of  the  estimate  of  xN  based  on  both 

the  correlation  between  all  of  the  variables  and  values  of 

the  previous  N-l  variables.  In  other  words,  -w  is  the 

estimate  of  x that  is  obtained  when  no  direct  knowledge  of 
N 

x is  available. 

N 

5.2  Linear  Spectrum  Extrapolation 

In  some  transform  coding  systems,  transmission 
bandwidth  is  saved  by  truncating  certain  spectral 
components,  notably  the  high  frequency  ones,  to  zero. 
(This  may  be  alternately  expressed  by  saying  that  these 
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components  are  transmitted  with  zero  bits.)  The  low 
frequency  components  are  retained  because  they  contain  more 
energy  than  the  high  frequency  components.  If  the  low  and 
high  frequency  components  are  correlated,  then  an  estimate 
of  the  high  frequency  components  can  be  obtained  by 
utilizing  this  correlation.  Since  higher  frequencies  are 
being  estimated,  the  estimation  technique  is  equivalent  to 
a spectrum  extrapolation.  If  the  transform  components  are 
also  samples  from  a gaussian  random  process,  they  can  be 
estimated  by  a Wiener  extrapolation  technique  developed  by 
Pratt  [3]  who  has  derived  the  following  linear 
extrapolation  operator 


T T -1 
W=CS  (SCS  ) 


(5.27) 


based  on  a minimum  mean-square  error  criterion,  in  which  C 
is  the  NxN  covariance  matrix  of  the  samples  and  the 
sampling  matrix  S has  the  form 


(5.28) 


where  is  the  (N-k)x(N-k)  identity  matrix,  and  0^ 
represents  k columns  of  zeros.  The  operator  provides  an 
estimate  of  k truncated  samples  based  on  N-k  known  samples. 


txrm 


Equation  5.26,  developed  in  the  last  section,  provides 
an  estimate,  yN , of  a single  vector  component  based  on  the 
known  samples  a^ ,a2 > . . • , 3^-1 • Note  that  this  also  is  a 
linear  estimation  of  a sample  from  a gaussian  random 
process.  The  estimate  obtained  by  this  method  is  identical 
to  that  obtained  from  eq.  5.^.7  for  the  special  case  k=l. 

To  demonstrate  this  equivalence,  consider  the 
following  partition  of  the  symmetric  covariance  matrix 


C= 


E \ 


— / 


(5.29) 


Note  that  D and  F are  also  symmetric  matrices.  For  the 
present  case,  let  D be  (N-l)x(N-l)  so  E is  Nxl  and  F is 
lxl.  Now  the  inverse  of  C can  be  found  in  terms  of  these 
sub-matrices  from  [4] 


/D' 1 + D" 1 E ( F - ET  D~ 1 E ) _]ET  D" 1 


E(F-ETD  1 


(5.30) 


\-  ( F~  ET£‘ 1 El 1 ET  D" 1 


(F-ETD_1 


In  the  notation  of  the  previous  section, 


T -1  -1 

rNN=<F-F2  E) 


(5.31) 
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Next  let  WN  be  an  (N-l)  vector  representing  the  Nth  row  of 
W,  which  was  previously  defined  in  eq.  5.27.  Then  the 
estimate,  y^,  using  this  vector  is 

yN  = WnTA  = (CST).)(SCST)~1A  = eV1^  (5.35) 

which  is  identical  to  eq.  5.33.  Therefore  eq.  5.26  can 

provide  a linear,  minimum  mean-square  error,  spectrum 
extrapolation . 


5.’  Covariance  of  Gaussian  Estimator 


The  conditional  covariance  of  the  estimate  found  in 
k e<?'  5,23  can  be  used  as  a measure  of  the  performance  of  the 


1 


estimate 


This  covariance  can  be  found  by  utilizing 


eq.  3.10  (without  the  subscript  m)  as  follows: 


cov  ( x 


f NttTP^)dxN 

J 

"ai 

• 

• 

f b 

• 

a..  , 

/ P(£> dxN 

N-l 

J aN 

yn 

Vl  yN  1 


Id  0 

0 0 


0 

0 


xNp(t)dx 


p(t)dxN 


N 


(5.36) 


Next,  denoting  the  (N,N)th  comoonent  of  this  matrix  as  eNN 
and  performing  he  integrations  reveals 


(bN+w) exp 

, , X T rNN  . . « 2 

- (e^  + w)  exD  — — (ajj+w) 

erf 

tTr‘v*>] 

+ 


1 

rNN 


- (y  +w) 2 


(5.37) 
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with  r and  w as  defined  in  eauations  5.16  and  5.24, 
NN 

respectively.  The  quantity  e represents  the  averaqe 

r i NN 

squared  error  that  results  when  x is  quantized  to  the 

N 

region  [a  ,b  ] and  an  estimate  of  xv,  is  obtained  based  on: 

N N N 

1.  the  quantization  region  [a  ,b  1 

N N 

2.  the  correlation  matrix,  C,  of  the  variables 

3.  the  previous  N-l  variables 

The  information  from  these  last  two  factors  is  contained  in 

w,  which  is  itself  a random  variable.  To  compare  e^,  to 

NN 

the  error  that  would  have  resulted  had  the  correlation 

information  been  unavailable,  it  is  necessary  to  remove  the 

conditioning  on  w.  This  can  be  done  by  averaging  e over 

NN 

all  values  of  w.  Because  w is  a weighted  sum  of  gaussian 

random  variables  (a  x,  to  a = x ),  w has  a gaussian 

11  N-l  N-l 

distribution.  The  parameters  of  this  distribution  are 

determined  next. 

The  mean  value  of  w can  be  found  from 


E { w} 


1 

rNN 


li'W 


0 


(5.38) 


The  variance  of  w is  then 


E { w2 } = E { (-i- 


N 

Er 


■NN  D 


= 1 N]  3 


2 

) } 


-^2  E{  <ITA)  ) 

'NN 


(5.39) 


where 
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E { w2  } = --2£TE{AAT}£  = -I2£TCnn £ (5.42) 

CNN  rNN 

where  C is  the  submatrix  of  C formed  by  deleting  the  Nth 
~NN 

row  and  the  Nth  column  of  C.  In  the  notation  of  Sec.  5.2, 
however,  C =D . Also, 


Then 


E { w2  } = ET  D 1E 


(5.45) 


so  that 

p(wj  =_____L___exp{"W^/(ETD  1 E)  ) (5.46) 

V2t(ETD-1  E) 

The  average  error  that  results  from  restoring  to  y^  a 
sample  that  has  been  quantized  to  the  interval  faN'bNl  is 
thus 


e 

Y 


p (w)  dw 


(5.47) 


Unfortunately,  a closed  form  solution  for  this  integral 


does 

not 

exist  when 

the  expression  for  has  the 

form 

given 

by 

eq . 5.37. 

Equation  5.47  therefore  must 

be 

evaluated 

numerically. 

The 

results  from 

one  such  numerical  evaluation 

are 

shown 

in 

table  5-1. 

The  integrations,  needed  to  calculate 

the  errors  listed  in  this  table,  have  been  performed  by 

means  of  a 16-point  gaussian-Hermite  quadrature  formula 

r 5 ] . The  quantization  regions  are  the  same  as  those 

calculated  by  Max  for  a unit  variance,  gaussian 
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TABLE  5-1 


NORMALIZED  MEAN-SQUARE  RESTORATION  ERROR  FOR 
A MAX-QUANTIZED,  GAUSS  I AN-MARKOV  PROCESS 
CORRELATION  = 0.9 
FILTERING  CASE 


Bits 

Quantization 

Interval 

Probability 

of 

Occurrence 

Normalized  Mean- 
Square  Restoration  Error 

Not  Using 
Correlation 

Using 

Correlation 

0 

i 

8 

+ 

8 

1 .00000 

1.000000 

0.190000 

1 

(0.000,+  oo  ) 

0 . 50000 

0.363344 

0.085441 

2 

(0.000,0.982) 

0.33685 

0.076896 

0.042723 

(0.982,+  oo  ) 

0.16315 

0.201243 

0.034435 

3 

(0.000,0.501) 

0.19165 

0.020687 

0.016412 

(0.501,1.050) 

0.16148 

0.024683 

0.017257 

(1.050,1.748) 

0.10662 

0.038160 

0.018193 

(1.748 ,+  oo  ) 

0.04023 

0.130557 

0.014966 

4 

(0.000,0.258) 

0.10188 

0.005543 

0.005140 

(0.258,0.522) 

0.09742 

0.005799 

0.005309 

(0.522,0.800) 

0.08871 

0.006373 

0.005667 

(0.800,1.099) 

0.07616 

0.007434 

0.006256 

(1.099,1.437) 

0.06047 

0.009386 

0.007122 

(1.437,1.844) 

0.04271 

0.013391 

0.008260 

(1.844,2.401) 

0.02445 

0.023947 

0.009205 

(2.401 ,+  oo  ) 

0.00818 

0.093356 

0.007747 

5 

(0.001,0.132) 

0.05251 

0.001451 

0.001420 

(0.132,0.265) 

0.05187 

0.001466 

0.001434 

(0.265,0.399) 

0.05067 

0.001502 

0.001466 

(0.399,0.536) 

0.04890 

0.001558 

0.001515 

(0.536,0.676) 

0.04653 

0.001636 

0.001584 

(0.676,0.821) 

0.04365 

0.001746 

0.001678 

(0.821,0.972) 

0.04026 

0.001893 

0.001795 

(0.972,1.130) 

0.03642 

0.002091 

0.001958 

(1.130,1.299) 

0.03222 

0.002367 

0.002175 

(1.299,1.481) 

0.02769 

0.002754 

0.002460 

(1.481,1.682) 

0.02295 

0.003325 

0.002841 

(1.682,1.908) 

0.01812 

0.004225 

0.003328 

(1.908,2.173) 

0.01332 

0.005764 

0.004036 

(2.173,2.504) 

0.00875 

0.008842 

0.004847 

(2.504,2.976) 

0.00467 

0.016968 

0.005638 

(2.976,+  oo  ) 

0.00146 

0.071316 

0.004866 
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distribution  [6].  For  a given  number  of  regions,  or, 


equivalently,  for  a given  number  of  bits,  this  choice  of 


t .gions  results  in  the  smallest  mean-square  quantization 


error.  This  error,  in  fact,  is  the  one  that  occurs  when 


the  restoration  corresponding  to  eq.  5.23  is  performed  and 


the  correlations  are  either  zero  or  simply  ignored. 


However,  when  successive  samples  of  the  distribution  are 


correlated  and  this  correlation  is  utilized,  the 


restoration  can  be  improved  and  the  mean-square  error 


decreased,  as  the  last  column  of  the  table  indicates.  The 


decrease  in  mean-square  error  is  seen  to  be  most 


significant  for  coarse  quantization  at  a low  number  of 


bits,  but  an  improvement  is  evident  for  every  quantizing 


region . 


The  specific  correlation  in  this  example  is  assumed  to 


be  due  to  a first-order  Markov  process  with  a correlation 


coefficient  of  0.9.  The  choice  of  a Markov  process  is  a 


convenient  one  because  the  restoration  of  a random  sample 


from  this  process  depends  only  on  the  sample  immediately 


preceding  <*>.  The  error  results  are  thus  independent  of 


the  number  of  samples,  N,  used  for  the  restoration,  and 


independent  of  the  size  of  the  correlation  matrix,  C.  When 


only  prior  samples  are  utilized  in  the  restoration  of  a 


<*>  By  the  definition  of  a Markov  process, 


P(XN  1 W--*'XN-1)=P(XN  1 Vl} 


■ ^ 


random  sample,  the  method  is  denoted  as  "filtering." 


Table  5-2  contains  similar  error  results  to  table  5-1, 
but  for  the  technique  of  "smoothing."  Smoothing  utilizes 
information  about  all  past  and  all  future  samples  to 
restore  a current  sample.  However  , in  the  case  of  a Markov 
random  process,  only  the  sample  immediately  preceding  and 
the  sample  immediately  following  the  current  sample  are 
needed  <**>.  The  use  of  this  added  information  from  future 
samples  significantly  reduces  the  resultant  mean-square 
restoration  error  for  any  choice  of  quantizing  region. 

The  results  of  tables  5-1  and  5-2  are  presented  in  a 
different  format  in  figures  5-1  and  5-2  to  illustrate  the 
total  mean-square  quantization  errors  that  occur  at  various 
quantizing  bit  assignments.  The  errors  are  calculated  bv 
using  the  information  from  the  tables  in  the  following 
equation 


M 


6=  52  e Pr{x  e (a  'bMl  } 
ym  N N N m 


(5.48) 


where  fa  ,b  ] is  one  of  the  M regions  to  which  x is 
N N m N 

quantized,  and  e is  the  error,  found  from  eq.  5.47,  for 

ym 


<**>  For  a Markov  process  it  is  also  true  that 


p(XN  | Xx  ,.  . . ,X^j  ,XN+1  ,.  . .)  PO^  | 
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NORMALIZED  MEAN-SQUARE  RESTORATION  ERROR  FOR 
A MAX-QUANTIZED,  GAUSSI AN-MARKOV  PROCESS 
CORRELATION  =0.9 
SMOOTHING  CASE 


Bits 

Quantization 

Interval 

Probability 

of 

Occurrence 

Normalized  Mean- 
Square  Restoration  Error 

Not  Using 
Correlation 

Using 

Correlation 

0 

( - OO  ,+  00  ) 

1.00000 

1.000000 

0 . 10  5000 

1 

(0.000  ,+  00  ) 

0 . 50000 

0.363344 

0.049106 

2 

(0.000,0.982) 

0.33685 

0.076896 

0.028041 

(0.982,  + oo  ) 

0.16315 

0.201243 

0.018555 

3 

(0.000,0.501) 

0.19165 

0.020687 

0.012428 

(0.501,1.050) 

0.16148 

0.024683 

0.012115 

(1.050,1.748) 

0.10662 

0.038160 

0.010920 

(1.748,+  oo  ) 

0.04023 

0.130557 

0.007810 

4 

(0.000,0.258) 

0.10188 

0.005543 

0.004501 

(0.258,0.522) 

0.09742 

0.005799 

0.004577 

(0.522,0.800) 

0.08871 

0.006373 

0.004720 

(0.800,1.099) 

0.07616 

0.007434 

0.004908 

(1.099,1.437) 

0.06047 

0.009386 

0.305144 

(1.437,1.844) 

0.04271 

0.013391 

0.005160 

(1.844,2.401) 

0.02445 

0.023947 

0.005431 

(2.401,+  OO  ) 

0.00818 

0.093356 

0.003924 

5 

(0.000,0.132) 

0.05251 

0.001451 

0.001356 

(0.132,0.265) 

3.05187 

0.001466 

0.001299 

(0.265,0.399) 

0.05067 

0.001502 

0.001333 

(0.399,0.536) 

0.04390 

0.001558 

0.001383 

(0.536,0.676) 

0.04653 

0.001636 

0.001444 

(0.676,0.821) 

0.04365 

0.001746 

0.001523 

(0.821,0.972) 

0.04026 

0.001893 

0.001638 

(0.972,1.130) 

0.03642 

0.002091 

0.001436 

(1.130,1.299) 

0.03222 

0.002367 

0.001719 

(1.299,1.481) 

0.02769 

0.002754 

0.001983 

(1.481,1.682) 

0.02295 

0.003325 

0.002267 

(1.682,1.908) 

0.01812 

0.004225 

0.002996 

(1.908,2.173) 

0.01332 

0.005764 

0.003345 

(2.173,2.504) 

0.00875 

0.008842 

0.003428 

(2.504,2.976) 

0.00467 

0.016968 

0.004066 

1 (2.976,+  OO  ) 

0.00146 

1 0.071316 

0.003735 
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I. 


SMOOTHING  CASE 


Fiqure  5-2.  Mean-square  restoration  error  for  a Max  quantized 
qaussian-Markov  process  with  correlation  factor  r. 


this  Quantization  region, 
significant  error  reductions 
utilizing  available  correlation 
fashion.  For  example,  when 
have  an  average  correlation  of 
one  bit,  utilizing  this 


The  errors  reveal  the 
that  can  be  obtained  by 
information  in  an  optimal 
samples  of  a random  process 
0.96  and  are  quantized  to 
correlation  in  an  optimum 


restoration  yields  a mean-square  error  eauivalent  to  that 
from  three  bit  Max  quantization  which  does  not  utilize  the 
correlation.  Figure  5-1  shows  the  filtering  case  and 
fig.  5-2  the  smoothing  case. 


5.4  Simulation  Results  for  Gaussian  Processes 

In  Sec.  5.1  a solution  is  described  for  acnieving  an 
optimal  nonlinear  restoration  of  quantized  gaussian  data. 
This  solution  has  now  been  applied  to  the  restoratioi^  of 
quantized  one-dimensional  random  signals  and 
two-dimensional  transform  domain  zonal-quant i zed  images. 
The  results  reveal  a decrease  in  mean-square  error  in  all 
cases.  However,  in  spite  of  the  error  reduction,  some 
images  exhibit  a degradation  in  subjective  quality  after 
being  restored.  Hence  a nonlinear  error  criterion  based  on 


the  human 

visual 

system  has 

been 

used 

in  place  of  the 

mean-square 

error 

function . 

Under 

this 

criterion  a 

subjective 

image 

improvement , 

as  well 

as  a 

numerical  error 

reduction,  are  obtained. 
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5.4.1  One-dimensional  Markov  Random  Process  Simulation 

To  demonstrate  the  utility  of  the  restoration 

procedure,  a randomly  generated  gaussian  Markov  signal  has 

been  quantized  and  restored.  A signal  having  unit  variance 

and  a Markov  correlation  of  0.95  is  first  generated.  A two 

bit  per  sample  Max  guantization  scheme  is  then  employed  to 

obtain  the  auantized  approximation  to  this  original  signal. 

The  quantized  signal  and  the  statistical  knowledge  about 

the  original  signal  are  used  as  inputs  to  the  nonlinear 

estimator  in  eq.  5.23.  The  restoration  is  performed  in 

block-lengths  of  16,  i.e.,  N=16  in  the  notation  of  this 

equation.  However,  the  inputs  to  the  equation  are  16 

quantization  intervals  rather  than  15  known  samples  and  one 

quantization  interval  as  the  estimator  requires.  Hence,  to 

satisfy  the  conditions  of  the  estimator,  a point  is  chosen 

within  each  of  the  15  intervals  as  an  initial  guess,  and  an 

estimate  of  the  sixteenth  sample  is  obtained  based  on  this 

guess.  This  procedure  is  then  repeated  for  each  quantized 

sample  in  turn,  with  the  guesses  successively  replaced  by 

their  calculated  estimates.  The  method  is  thus  recursive 

in  nature.  To  facilitate  the  convergence  of  these 

estimates  to  their  optimum  values,  the  initial  guesses  are 

chosen  to  be  the  Max  restoration  points  (when  the  sample 

correlation  is  ignored,  these  are  the  minimum  mean-square 

error  restoration  points).  It  has  been  experimentally 

determined  that  one  iteration  of  this  procedure  is 
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sufficient  for  convergence,  and  that  more  than  one 
iteration  provides  negligible  improvement. 

A typical  portion  of  the  results  from  the  above 
simulation  is  shown  in  fig.  5-3.  The  restoration  decreases 
the  mean-square  quantization  error  by  33%,  and  is  seen  to 
provide  a better  approximation  to  the  original  signal  than 
the  auantized  signal.  The  average  improvement  in 
mean-square  error  as  a function  of  quantizing  bit 
assignment  for  different  correlation  coefficients  is  shown 
in  fig.  5-4.  It  can  be  seen  from  this  graph  that,  as  the 
amount  of  correlation  in  the  Markov  process  approaches 
zero,  the  restoration  provides  no  error  improvement.  There 
is  no  improvement  as  the  number  of  quantizing  bits  becomes 
large  and  the  differences  between  the  original  signal  and 
the  quantized  signal  vanish.  There  is  also  r.o  improvement 
at  zero  bits  when  there  is  no  information  remaining  in  the 
quantized  signal  upon  which  to  base  a restoration.  Thus, 
the  above  procedure  represents  a viable  restoration 
technique  only  when  the  number  of  auantizing  bits  are  small 
and  the  input  samples  to  the  quantizer  are  correlated. 

5.4  2 Block  Transform  Zonal  Image  Coding 

The  conditions  which  were  placed  on  the  estimator 
derived  in  Sec.  5.1,  and  which  were  modified  experimentally 
in  the  preceding  paragraphs,  are  satisfied  by  the  zonal 
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transform  ceding  technique  for  images.  First,  transform 
samoles  typically  have  a gaussian  distribution:  each 
transform  sample  is  the  oum  of  a large  number  of  random 
variables  so  that  the  central  limit  theorem  can  be  invoked. 
Next,  for  all  transforms  except  the  Karhunen-Loeve , the 
transform  samples  are  correlated.  Finally,  to  achieve  a 
bandwidth  compression  or  a bit-rate  reduction,  some  of  the 
transform  samples  are  quantized  to  a small  number  of  bits. 
Hence,  because  all  of  the  necessary  conditions  for  a 
reconstruction  are  satisfied,  it  is  possible  and  practical 
to  restore  zonal  transform  coded  images. 


The  image  that  is  presented  to  a transform  coder  is 
assumed  to  be  in  the  form  of  a two-dimensional  array  of 
light  intensities.  (Each  point  of  this  array  is  known  as  a 
picture  element  or  "pixel”.)  If  the  image  is  presented  as 
a continuous  field  of  intensities,  however,  it  must  firs 
be  sampled  to  obtain  the  image  pixels.  A two-dimensional 
discrete  mathematical  transform  is  then  taken  of  these 
pixels.  The  transform  is  performed  over  the  entire  image 
or  over  subsections  of  the  image  known  as  blocks.  The 
transform  domain  samples  are  next  quantized  and  coded, 
either  for  storage  or  for  transmission  over  a channel  as 
fig.  5-5  shows.  At  the  receiver  the  samples  are  decoded, 
restored  to  reduce  the  quantization  effects,  and  inverse 
transformed  to  reconstruct  the  original  image.  The 
subsequent  paragraphs  discuss  these  operations  in  further 


r 


detail . 


5. 4. 2.1  Unitary  Transformations 


The  transform  operation  is  usually  performed  by  one  of 
a class  of  linear  unitary  operators  <*>  which  have  been 
recently  applied  to  image  coding  [7]  . A unitary  transform 
is  1 he  most  useful  because  it  preserves  the  "length"  or  the 
energy  of  a vector  component  during  the  transformation,  and 
because  it  transforms  every  orthonormal  basis  into  another 
orthonormal  basis.  To  interpret  this  concept  for  image 
codinq  purposes  reveals:  each  unitary  operator  can  be 

represented  Ly  a matrix  composed  of  orthonormal  basis 
functions  (generalized  spectral  functions);  each  component 
in  the  transform  domain  then  corresponds  to  the  amount  of 
energy  of  one  of  the  basis  (spectral)  functions  in  the 
original  (or  image)  domain.  Unitary  operators  transform 
sums  of  squares  into  sums  of  squares,  so  that  mean-square 
error  calculations  are  equivalent  in  either  domain.  This 
also  means  that  the  total  energy  in  the  original  and 
trarsform  domains  is  the  same.  The  unitary  transforms 
which  are  most  commonly  used  are 

1.  Karhunen-Loeve  transform 

2.  Cosine  transform 


<*>  A linear  operator  H such  that  HH*=H*H= I is  called  a 
unitary  operator,  where  * denotes  the  conjugate  transpose 
oT  the  complex  operator.  In  the  'eal  case,  H is  called 
orthogonal. 


' 


fc. 


1.  Fourier  transform 

4.  Slant  transform 

5.  Hadamard  transform 

6.  Laar  transform 
These  are  listed  in  order  of 
most  efficient  transform  to 
listed  in  reverse  order 

operation . 


efficiency,  i.e.,  from  the 
the  least.  They  are,  however, 
of  simplicity  and  speed  of 


The  Karhunen-Loeve  (K-L)  transform  results  in 

uncorrelated  transform  domain  components.  It  has  been 

found  that  quantizing  these  uncorrelated  components 

produces  a minimum  possible  mean-square  quantization  error 

[91.  Since  the  components  are  uncorrelated,  they  can  be 

quantized  and  restored  individually  without  loss  of 

performance,  so  that  the  restoratioM  method  of  Sec.  5.1 

provides  no  advantage  over  simple  Max  quantization  and 

restoration.  Less  efficient  transforms  result  in  greater 

quantization  error  than  the  K-L  transform  at  the  same  bit 

rate,  but  they  can  be  made  to  operate  much  faster.  Their 

transform  domain  components  have  some  statistical 

correlation  remaining  (if  the  components  in  the  image 

domain  were  ol  finally  correlated),  so  that  the  restoration 

method  of  Sec.  5.1  can  be  applied.  This  restoration 

reduces  the  mean-square  quantization  error.  The  remaining 

error  approaches  the  minimum  that  is  obtainable  by  means  of 

the  K-L  transform.  Thus,  almost  the  same  performance  as  a 
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F-L  transform  can  be  obtained  by  utilizing  a simpler 
transform  together  with  a more  complicated  reconstruction 
technique.  This  procedure  would  be  most  useful,  in 
practice,  when  an  image  coder  is  required  to  have  a simpler 
design  than  the  decoder.  For  the  image  coding  simulations 
and  reconstructions  in  this  chapter,  the  Hadamard  and  the 
Haar  transforms  are  utilized.  These  have  fast  and  easily 
implementable  algorithms,  and  also  provide  partially 
correlated  transform  domain  samples  which  can  then  be 
restored . 

5. 4. 2. 2 Zonal  Coding 

There  is,  in  general,  a uniform  distribution  of 
information  or  energy  throughout  an  image.  Unitary 
transforms,  although  preserving  the  same  total  amount  of 
energy  as  in  the  original  image,  rearrange  this  energy  and 
concentrate  it  in  a few  of  the  transform  domain  components 
[9] . This  concentration  is  achieved  because  the  transform 
makes  use  of  some  of  the  correlation  in  the  image.  Zonal 
coding,  which  takes  advantage  of  the  energy  concentration 
and  rearrangement,  entails  the  establishment  of  zones  of 
constant  bit  assignments  in  the  transform  domain 
corresponding  to  zones  of  approximately  constant  energy. 
The  bit  assignments  are  chosen  according  to  the  assumed 
variances — which  are  a measure  of  the  energy — of  the 
transform  domain  components.  A solution  for  this 
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assignment,  first  derived  by  Huang  and  Schultheiss  [10],  is 


where  V s the  matrix  of  transform  domain  variances  with 
components  j , b is  the  desired  average  number  of  bits  for 
each  component,  and  is  the  resultant  number  of  bits  for 
the  (i,j)th  component.  This  equation  provides  the  bit 
assignment  needed  to  quantize  gaussian  variables  based  on  a 
minimum  mean-square  error  criterion.  By  definition,  the 
number  of  bits  assigned  to  each  component  must  be  a 
nonnegative  integer.  Since  eq.  5.49  can  produce 
non-integer  and  even  negative  values,  tne  bit  assignment  it 
produces  must  be  adjusted  by  tr ial-and-er ror  techniques  to 
obtain  the  final  bit  assignment  for  each  component.  A 
typical  bit  assignment  for  a 16x16,  Haar  transformed,  image 
block  that  is  zonal  coded  with  an  average  of  one  bit  is 
shown  in  fig.  5-6. 


5. 4. 2. 3 Spatial  and  Transform  Domain  Correlation  Matrices 


To 

calculate  an  optimum 

bit 

assignment 

using  eq.  5 

.49, 

it  is 

necessary  to  know 

the 

var iance 

matfix  of 

the 

I 1. 

transform  domain  samples. 

This 

can  be 

derived  if 

the 

correlation  tensor  of  either  the  transform  domain  samples 


*■- 

or  the  image  domain  samples  is  known.  The  exact 


correlation  information  is  usually  unavailable,  but  the 
image  domain  correlation  can  often  be  accurately  modeled  as 
a Markov  process.  In  fact,  experimental  evidence  111,12] 
indicates  that  a reasonable  autocorrelation  function  for  a 
large  variety  of  p’.^torial  data  is  given  by 

C(x,x',y,y')=exp(-h  |x-x'|  -v  |y-y  |)  (5.50) 

(For  convenience  it  is  assumed  that  the  continuous  image, 
F(x,y),  has  zero  mean.)  This  function  can  be  used  to  model 
images  with  different  amounts  of  horizontal  and  vertical 
correlation  by  choosing  different  values  for  h and  v, 
respectively.  An  image  which  has  an  autocorrelation 
function  that  is  invariant  to  translation  is  said  to  be 
wide-sense  stationary.  The  autocorrelation  can  then  be 
rewritten  as 

C(Ax  ,Ay)  =exp ( -h  |Ax|  -v  |Ay| ) (5.51) 

This  autocorrelation  function  also  possesses  the  property 
of  horizontal  and  vertical  separability,  i.e.,  the 
correlation  between  any  two  points  of  the  image  is 
separable  into  the  product  of  horizontal  and  vertical 
correlation  functions. 

Now,  in  discrete  notation  for  an  NxN  sampled  image  F, 
the  correlation  between  any  two  pixels  can  be  represented 


by  a four-dimensional  tensor  as 


C(i, j,k,l)=E{F(i,k) F(j ,1) } (5.52) 

for  l£i,j,k,l<N,  It  is  sometimes  convenient  to  column-scan 
an  image  array  into  a data  vector,  f,  of  resultant  length 
N2  [13].  In  this  case  its  correlation  matrix,  , is 
N2  xN2  , and  can  be  represented  in  partitioned  form  as 


'-1,1  -1,2  ‘ ' * -1,N 

-2,1  -2,2  ‘ ‘ ‘ —2 , N 


(5.53) 


\ 


C C 
~N , 1 — N , 2 


. C, 


•N,N 


where  C is  the  correlation  matrix  of  the  (i)th  and  (j)th 
-if  3 

columns  of  F.  Under  the  assumption  of  wide-sense 
stationar ity , the  correlation  matrix  has  the  block  Toeplitz 
form 
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(5.54) 


& Sl-1 


Finally,  when  the  correlation  is  horizontally  and 
vertically  separable,  the  correlation  matrix  can  be  written 
in  direct  product  <*>  form  as 


(5.55) 


where  Cv  and  Ch  denote  the  NxN  correlation  matrices  of  the 
rows  and  columns  of  F,  respectively.  If  the  image  is  now 
considered  to  be  a sample  of  a Markov  process  with  a 
correlation  coefficient  of  v (0£v£l)  between  vertically 
adjacent  pixels  and  coefficient  h between  horizontally 
adjacent  pixels,  then 


<*>  The  symbol  8 denotes  the  matrix  direct  product.  The 
definition  employed  here  is  the  left  direct  product  [15]. 


•«  v < r 


where  sv  is  the  variance  of  the  pixels  in  each  column  and 
2 

is  the  variance  of  the  pixels  in  each  row.  After 
performing  a separable  unitary  transformation  H <*>  on  the 
image  F,  the  corresponding  row  and  column  correlation 
matrices  in  the  transform  domain  can  be  found  by 


<*>  For  a separable  transformation, 

H=H  B H, 

- -v  -h 

so  that 


C y=HGp  H 


,*T 


*/ 


= (i3v  8 SiHS,  ■ 5i><«v  8 3/T 

=(svsviT  > 8 <a,Si£T  > 


(5.58) 


The  matrices  and  can  now  be  used  to  find  the 

transform  domain  variance  matrix,  as 


V=V.V*T  (5.59) 

t-V- 


where  the  vectors  and  V consist  of 

elements  of  and  C^,  respectively. 


the  diagonal 


5. 4. 2. 4 Quantization  Levels 

Specification  of  the  variance  matrix  permits  the 
determination  of  an  optimum  bit  assignment.  However,  this 
provides  only  the  number  of  levels  to  which  each  component 
is  quantized  and  does  not  provide  the  locations  of  these 
levels.  For  the  analysis  presented  here,  the  placement  of 
these  levels  is  arbitrary,  but  is  chosen  to  be  the  Max 
quantization  levels  to  provide  the  minimum  mean-square 
quantization  error. 
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5. 4. 2. 5 Zonal  Coded  Images 


For  the  transform  image  coding  simulations,  the 
correlations  were  chosen  to  be  h=0.95  and  v=0.93.  These 
values  most  closely  modeled  the  correlations  of  the  images 
chosen  for  the  simulations,  and  resulted  in  the  smallest 
quantization  errors  and  reconstruction  errors.  Figure  5-7 
shows  the  three  "original"  images  which  were  used  for  the 
simulations:  each  image  is  an  array  of  256x256  pixels, 

with  the  intensity  at  each  pix»  1 quantized  to  one  of  256  (8 
bits)  grey  levels.  The  images  were  Haar  transformed  in 
blocks  of  16x16  pixels,  so  that  a 16x16  variance  matrix  was 
required  to  compute  the  bit  assignment.  Figure  5-8a 
displays  the  "girl"  image  after  being  zonal  Haar  transform 
coded  with  an  average  of  1.0  bits  according  to  the  bit 
assignment  shown  in  fig.  5-6.  (This  image  actually  shows 
the  Max  restoration  levels  after  the  image  has  been 
quantized  according  to  the  Max  decision  levels.)  Figure 
5-8b  shows  the  reconstructed  version  of  this  image  after  it 
was  restored  according  to  the  techniques  of  Sec.  5.1.  The 
restoration  technique  was  applied  recursively  to  the 
quantized  samples  since  it  is  only  capable  of  restoring  one 
sample  at  a time;  the  current  best  estimates  of  the  other 
samples  were  used  to  obtain  the  estimate  of  the  sample 
being  restored.  The  initial  estimates  of  the  samples  were 

chosen  to  be  the  Max  restoration  levels.  The  procedure,  in 
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Figure  5-7.  Origina 
simulations.  Each  im 
each  pixel  quantized 


MOON 


images  used  for  image  coding 
e consists  of  256x256  pixels 
8 bits. 


b.  Restored  image 


Figure  5-8.  Minimum  mean-square  error  restoration  of 
Haar  transformed,  one  bit  zonal  quantized  image. 


as  an  initial 


essence,  begins  with  the  Max  restored  image 
guess  and  attempts  to  improve  it  by  utilizing  the  known 
correlation  of  the  transform  domain  samples.  One  iteration 
of  this  restoration  procedure  was  utilized,  because  it  was 
experimentally  determined  that  more  iterations  were  of 
negligible  benefit.  The  mean-square  error  was  reduced  by 
10%  as  a result  of  the  restoration. 


The  "girl"  image  was  next  quantized  to  an  average  of 
0.5  bits  according  to  the  zonal  bit  assignment  shown  in 
fig.  5-9.  Figure  5-10a  shows  the  Max  restored  version  of 
this  coded  image  and  fig.  5-10b  shows  the  version 
reconstructed  according  to  the  technique  discussed  above. 
In  this  example,  the  reconstructed  version  has  19%  less 
mean-sauare  error.  The  "couple"  and  the  "moon"  pictures 
were  also  quantized  to  0.5  bits  and  restored.  The  results 
are  shown  i fig.  5-10.  The  resultant  mean-square  errors 
are  summarized  in  table  5-3.  The  restoration  technique  of 
Sec.  5.1  is  seen  to  provide  a significant  decrease  in 
mean-squavp  error  in  all  cases. 


5. 4. 2. 6 Visual  Coded  Images 

Subjectively,  the  reconstructed  images  of  figures 
5-8b,  5-10b,  5-10d,  and  5-10f  appear  to  be  much  smoother 
and  less  noisy  than  the  corresponding  Max  restored  images 
of  figures  5-8a,  5-10a,  5-10c,  and  5-10e.  However,  they 
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Figure  5-9.  Haar  transform  domain  zonal  quantizing  bit 
assignment  for  correlation  factors  h=0.95  and  v=0.93  and 
for  0.5  bits  per  pixel. 
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Quantized  image 


e.  Quantized  image 


f.  Restored  image 


Figure  5-10.  Minimum  mean-square  error  restoration  of 
Haar  transformed,  0.5  bit  zonal  quantized  image. 


TABLE  5-3 

NORMALIZED  MEAN-SQUARE  ERROR 
FOR  ZONAL  CODED  IMAGES 


Picture 

Error/Pixel 

after 

Quantization 

Error/Pixel 

after 

Restoration 

Improvement 

Girl 

1.0  bit/pixel 

1.024% 

0.922% 

9.9% 

Girl 

0.5  bit/pixel 

1.366% 

1.509% 

19.2% 

Couple 

0.5  bit/pixel 

1.793% 

1.412% 

21.2% 

Moon 

0.5  bit/pixel 

1.523% 

1.158% 

24.0% 

sometimes  appear  more  blurred,  as  is  very  evident  in 
comparing  figures  5-10e  and  5-10f.  An  improvement  in 
mean-square  error  apparently  does  not  correspond  to  an 
improvement  in  subjective  quality  in  all  cases.  Hence  an 
error  measure  is  required  in  which  numerical  results  match 
subjective  results. 

This  has  been  provided  by  modeling  the  error  measure 
after  the  human  visual  system.  Mannos  and  Sakrison  [14] 
have  derived  a nonlinear  error  criterion  which  achieves 
this  objective.  They  found  that  the  human  visual  system  is 
sensitive  to  approximately  the  cube  root  of  incident  light 
intensities.  It  is  also  most  sensitive  to  middle  spatial 
frequencies,  near  eight  cycles  per  degree  of  arc  subtended 
at  the  eye.  Hence,  to  apply  this  error  measure,  an  image 
has  been  processed  according  to  the  block  diagram  in 
fig.  5-11.  The  (i,j)th  component  of  the  filter  function 
shown  there  is  chosen  to  be 

Tij=(.05+.18525r)exp{- ( . 071 25r ) 1 ' 1 } (5.60) 

where 

r= ( i 2+ j 2) 1/2  (5.61) 

This  filter  was  applied  to  the  Hadamard  domain  seouencies 

of  the  image,  rather  than  the  Fourier  domain  frequencies, 
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11.  Coding  and  restoration  technique  for  a nonlinear  error  criterion 


because  of  the  close  similarities  between  seauencies  and 
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frequencies.  Figures  5-12a  and  5-12b  show  the  results  of 
the  procedure  for  a (1.5  bit  zonal  coded  Hadamard 
transformed  "moon"  image,  with  and  without  the 
reconstruction  operation,  respectively.  The  bit 
assignment,  shown  in  fig.  5-13,  was  calculated  by  choosing 
h=0.96  and  v=0.944  as  the  spatial,  cube-root  domain, 
correlations.  The  resultant  Hadamard  domain  variance 
matrix  wes  then  scalar  filtered  by  a multiplication  by  the 
square  of  the  filter  function  T.  The  mean-square  error  is 
reduced  by  5.5%  due  to  the  reconstruction  operation.  There 
is  also  a noticeable  subjective  improvement  after  the 
restoration.  Figures  5-12c,  5-12d,  5-12e,  and  5-12f 
display  similar  results  for  the  "girl"  and  "couole" 
pictures.  Table  5-4  summarizes  the  errors  for  the 
restorations  of  these  pictures. 

5. 4. 2. 7 Zonal  Coded  Color  Images 

The  restoration  technique  has  also  been  applied  in  an 
experiment  in  which  a color  image  is  encoded.  in  this 
experiment,  a color  image  is  first  transformed  to  the  YIQ 
color  coordinate  system  and  then  auantized  according  to  the 
bit  assignment  indicated  below  for  a typical  block  of  four 
pixels . 
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b.  Visual  restored  image 


d.  Visual  restored  image 


Figure  5-12.  Restoration  of  0.5  bit  Hadamard  transformed 
zonal  quantized  images  according  to  a visual  error 
criterion . 
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Figure  5-13.  A visual  quantizing  bit  assignment  in  the 
Hadamard  transform  domain  for  spatial,  cube-root  domain, 
correlation  factors  h=0.96  and  v=0.944,  and  for  0.5  bits 
per  pixel. 
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TABLE  5-4 

NORMALIZED  MEAN-SQUARE  ERROR 
FOR  VISUAL  CODED  IMAGES 


— 

Picture 

Er  roi/Pixel 
after 

Quantization 

Error/Pixel 

after 

Restoration 

Improvement 

Girl 

0.5  bit/pixel 

1.21% 

1.11% 

8.3% 

Couple 

0.5  bit/pixel 

1.94% 

1.79% 

7.5% 

Moon 

0.5  bit/pixel 

0.73% 

0.69%  1 

5.5% 
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Each  pixel  is  hence  coded  with  an  average  of  nine  bits  per 
pixel,  compared  to  an  original  coding  assignment  of  24  bits 
per  pixel.  The  restoration  technique  of  Sec.  5.1  provides 
a decrease  of  42%  in  mean-square  error  in  this  case  and  an 
improvement  in  subjective  quality. 


5.5  Summary 

Gaussian  data  which  are  correlated  and  which  have  been 
coarsely  quantized  are  amenable  to  being  restored  by  the 
techniques  outlined  in  this  chapter.  By  choosing  a 
suitable  error  criterion,  zonal  transform  coded  images  can 
be  analytically  and,  in  many  cases,  subjectively  improved 
so  that  they  more  faithfully  reproduce  the  details  of  an 
original  image.  The  restoration  techniques  have  been  found 
to  be  most  successful  when  the  transform  is  inefficient, 
i.e.,  when  some  statistical  correlation  remains  between  the 

transform  domain  samples,  and  when  the  quantization  is 
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coarse . 
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CHAPTER  6 

RESTORATION  OF  QUANTIZED  LAPLAC IAN  SAMPLES 


Laplacian  density  functions  arise  in  the  stochastic 
modelling  of  certain  kinds  of  communications  systems. 
Specifically,  the  output  of  a differential  pulse  code 
modulation  (DPCM)  system  can  often  be  modelled  as  a 
laplacian  random  process.  When  a laplacian  random  variable 
is  quantized,  as  in  a DPCM  system,  some  of  the  original 
information  about  the  variable  is  lost.  If  an  estimate  of 
the  continuous  laplacian  variable  is  then  made,  based  only 
upon  the  output  of  the  quantizer,  the  estimate  will  usually 
be  poor.  However,  if  the  characteristics  of  the  random 
process  are  known  and  are  utilized  also,  then  the  estimate 
can  be  improved  and  a reconstruction  of  the  original  signal 
will  be  attained.  The  estimation  equation  derived  in 
Chapter  3 (eq.  3.6)  provides  a means  for  accomplishing  this 
restoration.  Use  of  this  equation  requires  knowledge  of 


both  the  quantizer  structure  and  the  multidimensional 
probability  density  function  of  the  input  to  the  quantizer. 
An  approximation  to  a multidimensional  laplacian  density 
function  was  derived  in  Chapter  4.  This  density  function 
can  then  be  utilized  in  the  estimation  equation,  together 


with  arbitrary  quantization  parameters,  to  obtain  a 
restoration  of  quantized  laplacian  samples.  The  next 
section  Dresents  in  detail  the  solution  to  the  estimation 


ition  for  this  case.  The  results  of  that  solution  are 


then  applied  to  the  restoration  of  DPCM  coded  images, 


6.1  Laplacian  Quantization  Estimator 


A multidimensional  laplacian  probability  density 
function,  as  derived  in  Chapter  4,  can  be  approximated  by 


,-1/2 


N | N 


P(x) = — — exp{-V2£  Z<3ik\ 
2N/2  i=l  k = l 1K  K 


(6.1) 


where  C is  the  NxN  correlation  matrix  of  x,  gik  represents 
the  ( i , k) th  element  of  G as  given  by 


G=EA"1//2  ET 


(6.2) 


and  E and  A are  the  matrices  of  eigenvectors  and 
eigenvalues  of  C,  respectively.  This  density  function  can 


then  be  utilized  in  eg.  3.6  to  obtain  a minimum  mean-square 


error  estimate  of  a quantized  N-vector  of  laplacian 


samples,  x.  Now,  it  is  assumed  that  x is  quantized  to  a 


region  in  N-space,  R.  It  is  also  assumed  that  each 


component  of  x is  quantized  individually,  so  that  R is 


rectangular  (R  can  then  be  denoted  [a,b]).  Then  the 


estimate,  y,  of  x,  given  that  x G R,  is 


(6.3) 


N 

X,9ikxk 

k = l 


N 

£ 9ikxk 
k = l k 


} dx 


} dx 


A general  solution  to  these  integrals  exists,  but  is  rather 
complicated  to  state.  Special  cases  of  practical  interest 
are  considered  in  the  following  paragraphs. 


6.1.1  Scalar  Case 

When  the  components  of  x are  uncorrelated,  or  when 
they  are  restored  individually,  then  ea.  6.3  can  be 

decomposed  into  a product  of  one-dimensional  integrals 
which  can  be  solved  separately.  Considering  one  of  the 
components  of  x (for  a unit  variance), 


r b . 

/ 1xiexp{-V5'|  xi  | }dxi 

Llx 

f b. 

I 1exp{-V^| x^ | }dx^ 


( I b . | +^- ) exp{-\/I|  b,  | } - ( | a . | +— ) exp{-V2  I a . I } 

1 V2  1 1 s/2  1 

y _= 

1 b , a . 

— 1-  (exp{-\/2  I b . | ) -1)  - — --  (exp{-\/2  I a . I } -1) 

lb. | 1 la. | x 

i l 


(6.4) 
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If  a and  b,  are  chosen  to  be  Max  decision  levels,  then  yi 
i i 

is  the  corresponding  Max  restoration  level. 

6.1.2  Two-dimensional  Case 

In  two  dimensions,  eq.  6.3  yields  simultaneous 
estimates  of  two  quantized  and  correlated  laplacian 
variables.  The  solution  is  lengthy,  however,  so  the 
details  are  omitted  and  only  the  results  are  shown  here. 
To  simplify  the  notation,  let  the  exponential  terms  iri  the 
two-dimensional  version  of  eg.  6.3  be  written  as 

exp{- I fx^ -gx2  I - 1 gx^  -fx^  1} 

and  let  the  quantization  intervals  be  a<x^  <b  and  c<x2<d. 

The  constants  f and  g are  based  on  the  correlation  between 

x and  x . Also  let  s(*)  denote  the  sign  or  signum 
1 2 

function.  Then  the  numerator  of  y1  is 

{b [s (gd-fb) f +s ( fd-gb) g ] — 1 } ex  p { — I gd  — f b I — I fd-gb | } 

y ( num)  = — — —z  ’ ’ 

1 (s(gd-fb) f+s (fd-gb) g]2  [ s (gb-fd)  f+s  ( fb-gd)g] 

{b [s (gc-fb) f+s(fc-gb)g]-l}exp{-lgc-fbl-l fc-gb | } 

+ ' o 

[s  (gc-fb)  f + s ( fc-gb) g]2  [s(  fc-gb)  f + s (gc-fb) g] 

{a ( s (gd-fa) f+s ( fd-ga) g] -1} exp{- Igd-fa I -I fd-gal } 

[s (gd-fa) f+s ( fd-ga) g]2  [s( fd-ga) f+s (gd-fa) g] 
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{ a [ s (gc-fa) f+s  ( fc-ga;  g] -1} exp{- | gc-fa | - | fc-ga | } 


2 

(s (gc-fa) f+s ( fc-ga) g]  [ s ( fc-ga) f+s (gc-fa) g] 

[s  (gd-fb)  -s  (gd-fa)  ] gfl 2  [ |d  ( f2-g2  )/f  | +3 ] exp { - | d ( f2 -g2  ) / f I } 

( f2 -g2  ) 3 

[s  (gc-fa)  -s  (gc-fb)  ] gf2  [ I c ( f 2 -g2  ) / f I +3  ] exp{  - | c ( f2  -g2  )/f  I } 


[s  ( fd-gb)  -s  ( fd-ga)  ] fg2  [ |d  (f2  -g2  )/g  | +3]  exp{  - | d (f2  -g2  )/g  | } 


[s  (fc-ga)  -s  ( fc-gb)  ] fg2  [ I c ( f2  -g2  ) /g  | +3]  exp{-|  c(f2  -g2  )/g  | } 


[ s (gd-fb ) -s  (gc-fb)  ] g3  [ | b ( f*  -g2  )/g  | +1]  exp{-  | b ( f 2 -g2  )/g  | } 

(^-g2  )3 

[s(  fd-gb)  -s(  fc-gb)  J f3  [ I b ( f 2 -g2  ) /f  | +1  ] exp{-  | b ( t -g2  )/f  | } 

(f2  -g2  )3 

[s  (gd-fa)  -s  (gc-fa)  ] g3  [ | a ( f2  -g2  ) /g  | +1  ] exp{  - | a ( f2  -g2  )/g  | } 

( f2  -g2  )3 

[s  ( fd-ga ) -s(  fc-ga)  ] f 3 [ | a ( f2  -g2  )/f  | +l]exp{-  I a ( f2  -g2  ) /f  | } 

+ ( f2  -g2  )3 

(6.5) 


l 

Tne  denominator  is 
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ex  p { — | f b-qd | - | qb-fd | } 

y (denom)= 

1 [s (gd-fb) f+s ( fd-gb)q] [ s (gb-f d ) f +s ( f b-qd ) g ] 

exp{ - ! f b-gc | — | q b—  f c I } 

[s(qc-fb) f+s{fc-gb)q] [ s (gb-fc) f+s( fb-qc ) g] 
exp'- | fa-gd | -Iga-fd I } 

[ s (gd- fa) f +s ( fd-ga) g] [s(qa-fd)f+s( fa-gd ) g] 
exp{- | fa-gc I -Iga-fc I } 

4. — ‘ 

[s(qc-fa) f+s(fc-ga)g] [s(ga-fc)f+s(fa-gc)g] 

[s  ( fd-gb)  - s ( fc-gb)  ] f2  exp{  - | b (f2  -g2  )/f  I } b/f 
( f2  -g2  ) I b ( f 2 -g2  ) / f I 

[s  (gc-fb)  - s (gd-fb)  ] g2  exp{-  I b ( f 2 -g2  ) /g  I }b/g 
(f2  -g2)  I b ( f 2 -g2  ) / q I 

[s(fd-qa)  -s(fc-ga)  ] f2exp{-|a(f2  -g2  ) / f I }a/f 

+ — — — — — - — - — — — ■ — ■■  — - — — — 

( f2  -g2  ) I a ( f2  -g2  ) / f I 

[s(gc-fa)-s(gd-fa)  ] g2  exp{- | a ( f2  -g2  )/gl  }a/g 

+ — — — — — — — — — — — — — — — 

( f2  -g2  ) I a ( f2  -g2  ) /g  | 

[s  (gd-fb)  - s (gd-fa)  ] f2  exp{  - | d ( f;-  -g2  ) / f I }d/f 

+ 1 — — ■ — — — ■ ■■  1 - ■ ■ — — — — — — 

( f2  -g2  ) I d ( f 2 -g2  ) / f | 


[s(gc-fb)-s(gc-fa)  ] f2exp{-|c(f2  -g2  ) / f I } c/f 
( f2  -g2  ) I c ( f 2 -g2  ) / f I 


[s ( fd-gb)  -s ( fd-ga)  ] g2  expt- 1 d (f2  -g2  )/g  I } d/g 


(f2-g2)  |d(f2-g2)/gl 


[s  ( fc-gb)  -s(  fc-ga)  ] g^ exp{  - 1 c (f2  -g2  )/g  I } c/g 

+ — — (6.8) 

( f2  -g2  ) lc(f2-g2)/gl 


Hence , 


y,  (num) 


y ,=— 1 

1 y (denom) 


(6.7) 


The  other  component,  y^  , can  be  found  in  a similar  manner. 


Expressing  the  constants  f and  g in  terms  of  the 


correlation  between  x and  x„  reveals 

1 2 


VT+F  + Vl-t 


yj2(  1-r2) 


(6.8) 


Vl+r  - V1-l 


>j2(l-r2) 


(6.9) 


where 


r=E { x x } 
1 2 


(6.10) 


It  can  be  seen  that  as  the  correlation,  r,  approaches  zero, 


g also  approaches  zero  and  the  solution  (equations  6.5  to 


VT, ■**£.. 


6.7)  becomes  the  same  as  that  shown  in  eq.  6.4  for  the 
scalar  case.  For  non-zero  r,  equations  6.5  to  6.7  provide 
a minimum  mean-square  error  restoration  of  two  cor relaced 
laplacian  variables  which  have  been  quantised  to  the 
arbitrary  intervals  a£x^<b  and  c£x2.i^* 

6.2  Covariance  of  the  Laplacian  Estimator 

The  performance  of  the  estimator  found  in  equations 
6.5  to  6.7  can  be  analyzed  by  computing  its  conditional 
covariance  matrix.  Equation  3.10  contains  a general 
expression  for  finding  this  covariance  matrix  which,  in 
this  case,  must  be  solved  for  a two-dimensional  laplacian 
probability  density  function.  A general  solution  for  this 
covariance  matrix  has  been  obtained,  but  is  too  complicated 
to  be  shown  here.  Instead,  the  covariance  matrix  to  be 
derived  is  for  the  special  case  of  one-bit  quantization 
only. 


The  two  quantization  intervals  for  the  one-bit 
quantizer  are  chosen  to  be  [0,oo)  and  (-oo,0).  This  choice 
is  the  optimum  one  for  quantizing  individual  laplacian 
samples  according  to  a mean-square  error  criterion.  There 
are  now  four  possible  rectangular  regions,  R for 
m=l,2,3,4,  into  which  x=(x^,X2)  can  be  quantized.  Within 
each  region  there  is  a restoration  point  which  can  be  found 
by  means  of  equations  6.5  to  6.7.  The  restoration  points 
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are  dependent  on  the  correlation,  r,  between  and  x2  ; on 

whether  x^  is  quantized  to  a positive  or  a negative 

interval,  denoted  by  s(x^)$Q;  and  on  the  variance  of  x.  , 
2 

denoted  as  . Then  the  restoration  points  are 


1.  for  r s (x^ ) s (x2 ) >0 


<rl  (1+2  I r | ) s (x1 ) 

1 2 ( 2VT+T71  -Vl- 1 r | ) 


Y,  = 


(6.11a) 


y2  = 


°2  d+2  U I ) s (x2  ) 

2 ( 2vfi+| T -vT=TF|  ) 


(6.11b) 


2.  for  r s (x^ ) s (x2 ) <0 


/rm 

s(V 

/l-  I r | 


(6.12a) 


(6.12b) 


These  values  for  ^ can  now  be  substituted  into  eq.  3.10  and 
the  two-dimensional  integrals  evaluated  to  obtain  the 
covariance  matrix  for  each  of  the  two  cases  (assuming 

W111 


1.  for  rs (x^ ) s(x^ ) >0 
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cov{x  | x G R}= 


'Wl-r2-  (1-  I r | ) 2 I r | v4-r2+2  ( 1—  I r | ) 2 ' 

I 2v/T-r2-(l-|  rT)  4 [ 2 /l-r2  - ( 1-  | r | ) 1 

l |r|v/l-r2  + 2(l-|r|)2  2 v'l-r2  - ( 1- I r | ) 2 

\ 4 [ 2v£-72-  ( 1-  | r | ) ] 2/l-r2  - ( 1-  | r | ) / 


(1+2  I r |)  s(x,  ) 


( 1 + 2 I r | ) s (x,  ) 


( 1+2 | r I ) s (x  ) 


Vf  ( 2^1+1  r | -v/T-  I r | ) k/?  ( 2s/l+\r  I -^/l-  I r | ) \/2  ( 2/1+  I r | -/l-\  r 


( 1+2 | r | ) s ( x ) 
L 

\V2  ( 271+ 1 r 1 -yi-i  f n 


^ 4 ( | r I -2)  \/l-r2-2r2  + 9 
2(5+3  1 r |-4n/wT) 

(4-5|r |)/l-r2-8r2-2|r|-4 
v 4 (5+3  I r I - 4 v^T-  r 2 ) 


(4-5  | r | )\/l-r2-8r2-2  I r I -4^ 
4(5+31  r | -4%/l-r2  ) 

4(1  r|-2)\/l-r2-2r2+9 
2(5+3  I r | -4v/l-r 2 ) / 


(6.13) 


2 . for  rs (x  ) s (x  ) <0 
1 2 


The  probability  of  occurrence  of  each  quantization  region 
is  also  dependent  on  the  correlation.  These  probabilities 
are  found  by  integrating  the  two-dimensional  probability 
density  function  over  each  region.  The  results  are 


2.  for  r<0 


1 fr+r 


Pr  { s (x^  ) s (x2  ) >0 } =-^ 


(6.16a) 


411-r 


1 1 0 + r 


Pr  { s (x.  ) s (x2  ) < 0 } =-  - -J 

1 1 2 4%1-r 


(6.16b) 


Those  probabilities,  together  with  the  covariance  matrix, 


C,  for  the  original  vector  x, 


(6.17) 


and  the  expressions  for  the  restoration  points  listed  in 
equations  6.11  and  6.12,  can  be  used  in  eq . 3.13  to  obtain 


expression  for  the  total  mean-square  error  <*>: 


5=Tr{C  - L Y 4TPr(x  6 Rm>) 
m=l 


(3.13) 


which  becomes 


<*>  The  total  mean-square  error  is  the  same  whether  r is 
positive  or  negative,  so  only  the  results  for  r>0  are 
shown . 


uc^l; 


t-h  rn 


3r2+6r+(2-2r)Vl-r2  5r2+2r+2- (2-2r)Vl-r‘ 

8+8r-4Vl-r ^ 8 + 8 r - 4^1  - r ^ 


S=Tr(C  - 


5r2+2r+2-(2-2r)Vl-r2  3 r2+6r+ ( 2-2r )Vl-r7 

7 8+8r-4^1-r2  , 


S+8r 


-4\/l-r 


=Tr 


^8+2r-3r2+(6-2r)\/l-r‘ 

8+8r-4\/l-r2” 


3r2+6r-2+(2-2r 
8 + 8r-4\/l-r^ 


\ 


3r2+6r-2+(2-2r  )Vl-r2 
-4Vl-r2 


2 ■ " --'Vl-r2 


8+8r- 


8+2r-3r  + (6-2r) 
8+8r-4\/l-r2 


8+2r-3r2-(6-2r)Vl-r2 
4+4r-2  


(6.18) 


■ 

This  error  is  plotted  in  fig.  6-1  as  a function  of  the 
correlation  r.  Also  shown  is  the  mean-square  error  that 

results  when  the  correlation  between  any  two  quantized 
laplacian  samples  is  not  utilized  in  the  reconstruction  of 
these  samples.  The  use  of  the  correlation  information  is 
seen  to  substantially  reduce  the  resultant  mean-square 
error.  This  holds  true  independent  of  the  sign  of  the 
correlation,  but  is  only  true  when  there  is  at  least  a 
moderate  amount  of  correlation.  The  next  section  shows 
that  deltamodulation  image  coders  satisfy  this  requirement, 
so  that  a reconstruction  of  the  coded  image  can  then  be 

I 
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obtained . 


6.3  DPCM  and  Deltamodu] ation  Image  Coding 

Differential  pulse  code  modulation  (DPCM)  is  a common 
technique  for  achieving  bandwidth  compression  in  digital 
systems.  In  this  method  differences  between  successive 
signal  samples,  rather  than  the  signals  themselves,  are 
transmitted.  Compression  occurs  because  adjacent  samples 
are  often  very  similar,  and  transmitting  only  signal 
differences  removes  some  of  this  redundancy.  A block 
diagram  of  a typical  DPCM  system  is  shown  in  fig.  6-2.  An 
essential  component  of  this  system  is  the  quantizer, 
because  it  permits  a bandwidth  compression.  The  coarser 
the  quantization,  the  greater  the  compression,  but  also  the 
greater  the  degradation  of  the  reconstructed  signal.  To 
minimize  this  degradation,  a reconstruction  must  utilize 
all  of  the  knowledge  that  is  available  about  the  signal, 
such  as  the  quantization  levels,  the  signal  distribution, 
and  any  correlation  which  remains  after  the  differencing 
operation.  Section  6.1  presented  a means  for  achieving  an 
optimal  solution  to  this  restoration  problem.  This  section 
applies  that  solution  to  the  reconstruction  of  DPCM  coded 
images . 

In  the  DPCM  system  of  fig.  6-2,  the  difference  signal 
d(k)  is  formed  by  subtracting  a prediction  of  the  current 
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signal  from  the  current  signal  itself.  The  first-order 
predictor  shown  in  the  figure  can  be  characterized  by 


specifying  the  constants  a^  and  a^  . These  constants  are 
usually  chosen  to  minimize  the  mean-square  prediction 
error.  The  prediction  error  can  be  written  as 


£=E{  [x(k)-(alX(k-l)+a0)]  } 


(6.19) 


This  expression  can  be  minimized  by  simultaneously  solving 
the  two  equations 


-2m+2a1m=0 


(6.20a) 


and 


1 : 


r — =2an  s -2 r sz  + 2a  m=0 
dax  1 0 


(6.20b) 


where 


m=E { x ( k ) } 
s2  = E { [ x ( k) ] 2 } 


(6.21) 


(6.22) 


and 


E{x(k) x(k-l) } 


(6.23) 


The  results  are 


ms2 (1-r ) 


a°”  *2  m2 

s -m 


(6.24) 


2 2 
rs-m 

a,  = — — — 

1 o2  _2 

s -m 


(6.25) 


The  reconstruction  unit  in  the  block  diagram  of  the 
receiver  in  fig.  6-2  is  a device  which  attempts  to  reduce 
the  effects  of  the  quantizer.  The  particular  form  of  this 
device  is  based  on  a priori  knowledge  of  the  quantizer  and 
the  statistics  of  the  quantizer  input.  This  input  is  the 
difference  signal,  which  has  been  found  to  have  a laplacian 
distribution  [1] , described  by 


P ( x ) — — ■* — exp(-|  x I / 2 cr  ) 
y/2  <j 


(6.26) 


The  distribution  of  the  difference  signal  for  the  "girl" 
image  (see  fig.  6-6)  is  shown  in  fig.  6-3.  It  can  be  seen 
that  a laplacian  could  model  this  distribution  quite  well. 
However,  it  has  been  found  that  the  difference  signals  of 
images  are  correlated,  so  that  an  accurate  statistical 
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representation  of  them  must  also  account  for  this  fact. 
Figure  6-4  shows  the  actual  two-dimensional  distribution  of 
the  DPCM  coded  "girl".  For  tnis  image,  the  average 
correlation  for  adjacent  difference  samples  has  been 
measured  as  0.4.  These  samples  car.  then  be  modeled  by  a 
correlated  two-dimensional  laplacian  density,  written  as 


P(x,y)= 


2<t  a \l  1-r 
x y* 


5 exp{- 


\l  2(1-0 


ax_b 

4. 

ay  bx  | 

a a 
x y 

T 

a a 
y X 1 

)}  (6.27) 


where 


*x2=E{x2} 

(6.28) 

ay2=E{y2} 

(6.29) 

E{xy} 

r = 

a 

(6.30) 

x y 

a=^l+r  + yjl-r 

(6.31) 

h=yjl+r  - yjl-r 

(6.32) 

Figure  6-5  contains  a plot  of  this  density  function  for 

r=0 . 4 and  o = a . 

x y 


The  two-dimensional  laplacian  distribution  is  seen  to 
provide  an  accurate  model  for  the  DPCM  samples.  These 
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samples  are  quantized  before  thev  are  transmitted  throuqh 
the  channel.  A minimum  mean-square  error  restoration  of 
the  quantization  can  be  obtained  by  utilizinq  equations  6.5 
to  6.7.  This  solution  is  complicated,  but  can  be 
simplified  for  the  case  of  deltamodulation  (one  bit 
quantization).  In  this  case  x and  y are  each  auantized  to 
the  intervals  [0,°O)  or  (-00, 0),  or  equivalently,  as 
positive  or  neqative.  The  results  are  contained  in 
equations  6.11  and  6.12. 

Applyinq  this  restoration  to  the  quantized  imaqe  in 
fiq.  6-6a  results  in  a mean-square  error  reduction  of  12%. 
Subjectively,  as  the  restored  image  in  fiq.  6-6b  shows, 
there  is  less  apparent  noise  and  more  discernible  detail. 
Equations  6.5  to  6.7  have  also  been  utilized  directly  to 
restore  samples  quantized  to  general  regions. 
Specifically,  the  results  were  applied  to  the  two  and  three 
bit  coded  images  shown  in  figures  6-6c  and  6-6e, 
respectively.  The  reconstructed  images  in  figures  6-6d  and 
6-6 f exhibit  bo  :h  a reduction  in  mean-square  error  and  a 
subjective  improvement  in  quality.  The  subjective 
improvement  is  less  apparent  in  these  pictures,  however, 
because  the  quantization  itself  is  less  noticeable.  Table 
6-1  shows  the  quantization  intervals  that  were  used  to  code 
these  images  [21 . The  intervals  were  chosen  to  minimize 
the  mean-souare  quantization  error  that  occurs  when  an 
individual  laplacian  sample  is  quantized. 
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Figure  6-6.  Minimum  mean-square  error  restoration  of  DPCM 
encoded  images  using  two  adjacent  pixels. 


TABLE  6-1 


QUANTIZATION  INTERVALS  FOR  SIGNALS  WITH  A LAPLACIAN 
DISTRIBUTION  CHOSEN  ACCORDING  TO  A MINIMUM 
MEAN-SQUARE  ERROR  CRITERION 


Bits 

Quantization 

interval 

Restoration  Point 

0 

- oO 

0.0 

1 

0.0 

oO 

0.707 

2 

0.0 

1.102 

0.395 

1.102 

OO 

1.810 

3 

0.0 

0.504 

0.222 

0.504 

1.181 

0.785 

1.181 

2.285 

1.576 

2.285 

OO 

2.994 

4 

0.0 

0.266 

0.126 

0.266 

0.566 

0.407 

0.566 

0.910 

0.726 

0.910 

1.317 

1.095 

1.317 

1 .821 

1.540 

1.821 

2.499 

2.103 

2.499 

3.605 

2.895 

3.605 

OO 

4.316 

the 
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The  image  reconstructions  shown  in  fiq.  6-6  applied 
techniques  summarized  in  equations  6.5  to  6.7  by  usinq 
the  quantized  information  about  a previous  sample  to 
restore  a current  sample.  By  solving  eq.  6.3  in  three 
dimensions  rather  than  in  just  two,  quantized  information 
about  two  other  samples  can  be  used  in  the  reconstruction 
of  a third  quantized  sample.  The  DPCM  image  samples  are 
thus  modelled  by  a correlated  three-dimensional  laplacian 
density  function.  A form  of  this  three-dimensional 
solution  was  then  employed  to  obtain  the  reconstructed 
image  shown  in  fig.  6-7a.  In  this  imaqe,  the  sample 
immediately  preceding,  and  the  one  immediately  following, 
the  sample  being  restored  were  used  in  the  reconstruction. 
The  resultant  image  has  39%  lower  mean-square  error  than 
the  one-bit  quantized  image  shown  in  fig.  6-6a.  Figure 
6-7b  presents  a similar  result,  except  the  samples 
immediately  above  and  below  the  current  sample  were 

utilized  in  the  reconstruction.  In  this  case,  the 
mean-square  error  improvement  is  18%.  Both  image 

reconstructions  also  reveal  a distinct  subjective 
improvement — there  is  a reduction  in  visual  noise  in 
constant  luminance  areas  of  the  image,  together  with  a 


t/  . 


decrease  i 

n slope 

overload  <*>  at  edges  with 

in  the  image. 

Thus 

the  technique  described 

above 

provides  an 

effective 

method 

for  restoring 

DPCM 

coded  images, 

particular 

ly  when 

the  quantization  is 

coar  se 

and  the  DPCM 
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Restoration  using  three  horizontal  pixels 


Figure  6-7.  Minimum  mean  -square  error  restorations  of 
one  bit  (deltamodulation)  encoded  image  using  three 
adjacent  pixels. 


samples  are  correlated. 
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CHAPTER  7 

QUANTIZATION  AND  RESTORATION  OF  RAYLEIGH  SAMPLES 


The  last  two  chapters  have  shown  that  the  nonlinear 
Quantization  restoration  equations  of  Chapter  3 can  be 
solved  for  both  gaussian  and  laplacian  quantizer  inputs. 
This  chapter  extends  these  results  to  quantizer  inputs 
which  can  be  described  by  Rayleigh  probability  density 
functions.  It  will  be  shown  that  a Rayleigh  distribution 
can  accurately  model  the  intensity  distribution  of  an 
image,  so  that  quantized  images  can  be  reconstructed  by  the 
techniques  developed  herein. 

7.1  Rayleigh  Densities  in  PCM  Image  Coding 

In  the  PCM  (pulse  code  modulation)  coding  of  an  image, 
the  image  is  first  sampled  at  an  array  of  points  known  as 
pixels.  The  value  assigned  to  each  pixel  is  a measure  of 
the  light  intensity  at  that  point.  Each  light  intensity  is 
then  quantized  and  coded,  for  either  digital  transmission 
or  storage.  Because  it  is  proportional  to  the  square  of 
the  magnitude  of  an  electric  field,  light  intensity  is  a 
real  and  non-negative  quantity.  The  distribution  of  light 
intensities  for  an  image  can  be  characterized  by  a 
histogram  ranging  from  black  (zero)  to  white.  For  most 
ratural  images,  there  are  many  more  dark  pixels  than  bright 


t 


pixels,  and  the  envelope  of  the  histogram  tends  to  fall  off 
exponentially  at  higher  brightness  levels  [1] • Image 
histograms  measured  by  Kretzmer  [2]  and  Stockham  [3] 
exhibit  these  properties.  Figure  7-1  shows  that  the 
histogram  of  the  "girl"  image,  shown  originally  in 
fig.  5-7,  also  possesses  these  characteristics. 

In  the  past,  image  intensity  distributions  have  been 
well-modelled  by  Rayleigh,  log-normal,  and  exponential 
probability  density  functions  [4).  In  this  chapter,  a 
variation  of  a Rayleigh  distribution  extended  to  many 
dimensions  is  utilized  to  model  image  intensities.  A need 
for  a multidimensional  distribution  arises  because  of  the 
usefulness  of  including  pixel-to-pixel  correlations  in  the 
model.  A correlated  multidimensional  Rayleigh 
distribution,  derived  in  Chapter  4 as  eq.  4.34,  can  be 
rewritten  as 


P(x)  = l£l_1/2TT(-i-)U(-i-)exP{“3-T£  ^ (7.1) 

i=l 


where  C is  the  correlation  matrix  of  x and  the  vectors  h^ 
are  a function  of  the  eigenvalues,  X- , and  eigenvectors, 

I ‘ 

e.,  of  C.  If  the  eigenvectors  are  arranged  column-wise 
into  a matrix  E having  elements  e.^,  then  the  components  of 

[ *, 

the  vectors  h^  can  be  found  from 
fc. 


t cn 


(7.2) 


ik 


N e.  . e,  . 
= E_imi 


FiV\ 
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The  average  correlation  between  adjacent  pixels  for  an 
image  such  as  the  "girl"  has  been  measured  to  be  0.95.  A 
one-dimensional  marginal  probability  distribution  of 
eq.  7.1  for  N=2  and  a correlation  of  0.95  is  shown  in 
fig.  7-2  (an  expression  for  this  marginal  distribution  can 
be  found  in  eq.  4.45).  This  distribution  is  seen  to 
provide  an  accurate  model  for  the  intensity  distribution  of 
fig.  7-1.  A two-dimensional  histogram  of  the  "girl"  image, 
obtained  by  plotting  pairs  of  adjacent  pixel  intensities, 
is  shown  in  fig.  7-3.  This  histogram  can  be  closely 
modelled  by  the  correlated  two-dimensional  Rayleigh  density 
function  shown  in  fig.  7-4.  (Only  the  positive  quadrants 
of  figures  7-3  and  7-4  are  shown  since  both  the  histogram 
and  its  Rayleigh  model  are  identically  zero  elsewhere.) 
This  statistical  model  for  image  intensities  can  then  be 
used  for  the  restoration  of  quantized  and  coded  values  of 
these  intensities. 


7.2  Estimation  of  Quantized  Rayleigh  Samples 

In  order  to  obtain  a restoration  of  quantized  Rayleiqh 
samples,  such  as  those  that  result  from  the  PCM  coding  of 
an  image,  eq.  3.6  must  be  solved  using  eq.  7.1  as  a 

statistical  model  of  the  underlying  Rayleigh  random 
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process.  Unfortunately,  an  explicit  solution  to  eq.  3.6 
does  not  exist  for  this  choice  of  input  distribution,  p(x) . 
However,  solutions  can  be  found  for  several  special  cases 
of  quantization  regions  and  correlation  matrices.  Also,  an 
approximate  solution  can  be  obtained,  using  a method 
similar  to  Curry's  [5],  which  is  valid  for  very  fine 
quantization . 

The  approximate  solution  can  be  found  by  means  of  a 
power-series  expansion  of  p(x).  By  retaining  terms  of  the 
expansion  up  to  fourth-order,  the  approximation  is  accurate 
enough  to  restore  the  outputs  of  a wide  variety  of 
quantizers.  A sufficient  condition  for  the  approximation 
to  hold  is  that  the  size  of  the  quantization  intervals  must 
be  small  compared  to  the  variance  of  a component  quantized 
to  one  of  these  intervals.  A multidimensional  Taylor 
series  is  used  to  expand  p(x)  about  the  midpoint,  cj_,  of  a 
rectangular  quantization  region,  R.  Let  the  region  be 
defined  as 


R={x  | x e [ a , b] } (7.3) 

where  a,  b,  and  x are  N-vectors,  so  that  its  midpoint  is 

g=(a+b)/2  (7.4) 

The  size  of  the  quantization  region  can  be  expressed  in 
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terms  of  a vector  of  half-widths  as 


d= (b-a) /2 


(7.5) 


Taylor-ser ies  expansion  of  p(x)  about  $ is  then 


&p(*>  t Tap(*)  # . 

P(X)=P(°I)+  -7—  ;x-2)+2  (-"2)  ^xT  -~2 


N N N d p(x) 

4L  £ 'VV'VV'W*  ’ ' 

6 i=l  j=l  k=ldx. dx. dx 

1 3 K a 


. (7.6) 


where  p(x)  is  given  by  eq.  7.1.  The  terms  of  this 
expansion  are  evaluated  individually  as  follows: 


?,P(x)  N 1 Mh.x)  -l  | 

=p ( x ) y.  h,  — — + — — ~ £ * 

Sx  “ i = l_1  hTx  U(hTx) 


(7.7) 


where  &(•)  represents  the  Dirac  delta  function,  defined  as 
the  derivative  of  tie  unit  step  function  U(*)»  and 


£=lhil  hi2  • • • hiN> 


(7.8) 


This  vector  represents  the  (i)-h  row  of  the  matrix  H. 
Next,  define  the  vector  w such  that  its  components  are 


1 


so  that  eq.  7.7  can  be  rewritten  as 


dp(x) 

dx 


= P(g)  (HTW-C  1q) 

a 


The  second  derivative  can  be  computed  in  a similar 
as 


iT, 


dx 


C*P(X) 

dx 


=P(g)  [ (HTw-C-1g)  ( HT  w - C”  1 q ) T - C" 1 - fT  KH  ] 


where  W is  a diaqonal  matrix  havinq  components 

1 8'  ( hTg ) 5(hJ’a)U(h^q) 

(h^g)  U(h^g)  [U(h^)]2 

Now  let  z=x-c[.  Then  the  denominator  of 
restoration  equation  (eq.  3.6),  after  retaininq  the 
three  terms  of  the  expansion,  becomes 


(7.10) 

manner 


(7.11) 


(7.12) 

the 

first 


/b  f d T -1  T 

p(x)dx  p(£){l+(H  w-C  <£)  z 


+2ZT[ (HTw-C_1g) (HTw-C_1g)T-C_1-HTWH] z}dz 


1 « . di 


■ p(g)  V (1  + 2 £ Bii-T 
1 = 1 


(7.13) 


where 


rn  _1  T -1  T “I  T 
B=((H  w-C  a) (H  w-C  g)  ~C  -H  WH] 


(7.14) 


and 


N 

V=TT  (2d.) 

i=l 


(7.15) 


The  quantity  V is  the  volume  of  the  rectangular 
quantization  region  R,  and  B is  a matrix  based  on 
parameters  of  this  quantization  region. 

Next,  the  numerator  of  ea.  3.6  can  be  approximated  by 
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/ xp(x) dx=  / (g+z)  p(g  + z)  dz 

J a J -d 

f—  f — T -1  T 1 T 

~g  / p(x)  dx  + / zp(£)  { 1+ (H  w-c  q)  z+jz  Bz}dz 


N d 


P<2)Vtl+j.E  Bii~  + 7(HTw-C_1£)D] 


i = l 


3 3 


(7.16) 


where  the  diagonal  matrix  D has  components 


D . 
.1 


2 


(7.17) 


Thus,  the  conditional  mean  estimate  of  a vector  x, 
quantized  to  the  region  R,  is 


(HT  w-C_1g)D 

y=a  + — \ 

3 + jTr{BD} 


(7.18) 


This  estimate  is  seen  to  consist  of  a correction  term  added 
to  the  midpoint  of  the  quantization  region.  Since  the 
estimate  is  a function  of  a particular  quantization  region, 
the  estimation  process  is  a nonlinear  operation.  The 
estimation  results  are  valid  for 


2d.  <s. 
i i 


(7.19) 
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for  i=l,2,...,N,  where  s^  is  the  variance  of  the  (i)th 
component  of  x.  It  has  been  experimentally  determined  that 
an  image  must  be  coded  with  more  than  four  bits  per  pixel 
for  the  solution  given  by  eq.  7.18  to  be  useful. 

For  very  coarse  quantization,  the  relation  in  eq.  7.19 
is  not  satisfied,  so  another  solution  to  eq.  3.6  must  be 
obtained  which  is  valid  f this  situation.  A solution  can 
be  found  by  utilizing  a metnod  developed  in  Chapter  5.  The 
essence  of  this  method  is  that  the  components  of  a 
quantized  vector  are  restored  individually,  based  on 

1.  the  quantization  interval  of  the  component  being 
restored ; 

2.  estimates  of  the  remaining  quantized  components 
(not  the  quantization  intervals  of  the  remaining 
components) . 

The  initial  estimates  can  be  chosen  to  be  the  midpoints  of 
the  quantization  intervals.  Alternately,  they  can  be 
chosen  as  the  solution  to  the  one-dimensional  version  of 
eq.  3.6,  which  can  be  solved  exactly.  This  latter  choice 
is  the  more  accurate  one,  and  is  the  choice  that  will  be 
made  here.  As  new  estimates  of  the  components  are  found 
using  the  method  of  Chapter  5,  they  are  used  to  replace  the 
initial  estimates.  The  technique  is  thus  iterative. 

Keeping  these  concepts  in  mind,  consider  a typical 
component  of  eq.  3.6.  This  component  can  be  solved 
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where  s£  is  the  variance  of  xR  . Unfortunately,  eq.  7.20 
cannot  be  integrated  so  easily.  When  the  expression  for 
p(x),  as  defined  by  eq.  7.1,  is  substituted  into  this 
equation,  the  estimate  becomes 


/b  • N m m 1 T 1 

xiTT-i-U  (-i-}  exp*~2-  £ t}dxi 

/b*  N m m IT  “1 

JT-i-U(-i-)exp{"Z-  - -}dxi 

a.  k = l 


(7.23) 


where 


tT=[t1  t2  . . . ti_1  Xi  ti+1  . . . tN]  (7.24) 


Proceeding  as  in  equations  5.16  to  5.19,  common  factors  can 
be  cancelled  from  both  the  numerator  and  denominator, 
leaving 


fb.  N -1  2 

/ lxiJT(hkixi+Hk)u(hkixi+Hk»exPf; — (riixi+Bi'  )dxi 
Ja.  k = l 2r  jj 

i_  fb . N -1  2 

1 IT (h  x.+Hk)U(hkixi+Hk)exp{ (riixi+Ri>  *dxi 

J a . k = 1 2r  ■ • 


where  r _ is  the  (i,j)th  element  of  C , 


r 


N 


^ fcj 


and 


R.  = E r.  .t. 

1 j = l ^ 3 


j*i 


Making  the  substitution 


z=- 


V^Ti 


(r.  . x.  +R.  ) 
111  1 


yields 


/v  _ N 

(J27.  . z-R.  )JT“ 
u^  11  1 k-l 


y.  =■ 

i 


\i 


r.  . 
L 11 


S Phlz-Ri'+S 


exp{-zz }dz 


li 


Tv  N 

h . 

/ JT 

/ u k = l 

— (\l 2r  . . z-R.  ) +H, 
r » n i k 

ii 

exp{ 

where 


max 
u=  k 


h .>0 
ki 


V2 


1 R.  H,  frTT 

( r . . a . +R. } , i S-|-ii 

T7  11  1 1 \/2c77  h,.V  2 


1 1 


ii  ki’ 


mm 
v=  k 

h <0 
ki 


1 R.  H,  frTT 

— (r..b.+R.),  ~ - -JS-J-li 

r,,  11  1 1 V2TT  h.  ,V  2 


11 


ii  ki 


(7.26) 


(7.27) 


(7.28) 


(7.29) 


(7.30a) 


(7.30b) 


The  expressions  which  multiply  the  exponential  terms  in 
eq.  7.29  are  seen  to  be  Nth-order  polynomials.  Writing 


"T’ * ±r. 


'■**!  ,‘V -XT  r 


id. 


174 


them  as  such,  reveals  that 


_ fv  N 

I z£  g,  zkexp{-z2}dz 

J u k = 0 k 


P. 


y*- 


1 /*  v N 

vn  / . 5 


(7.31) 


r.  . 


E g z exp{-z  }dz 


li 


The  integrations  can  finally  be  performed  to  yield  <*,**> 


num 


. =yj2 


|N±1| 

yjieitw 

k = l 2k 


N 


k-2j 


-,v 


-exp{-z2 } X!  g.  k!  ! 2^ 

i e = a K -;  = a 


k=0  K j=0  23  ( k-2 j ) ! ! 


(7.32a) 


denom.=\/r . . 

* li 


N I 


k=0 


N 


M 


k-2 j-l 


iV 


-exp{  -z^  } £ 9.  (k-D  ! ! E 


i—4 

k = l K 


j = 0 P.3  ( k-2  j-l)  ! 1 


(7.32b) 


so  that 


<*>  The  notation  |_xJ  denotes  the  smallest  integer  < x. 
<**>  The  notatior.  k!l  denotes  the  product 

!k  * ( k — 2 ) - • *3*1  for  k odd 

k*(k-2)*  • *4*2  for  k even 


num. 


(7.33) 


R. 

1 denom.  r. . 

1 i 

When  C=_I  or,  equivalently,  when  x^  is  being  restored 
independently  of  the  other  components  of  x,  the  estimate 
becomes  yi=ti,  as  given  by  eg.  7.22.  When  C=J[,  the 
estimate  in  eq.  7.33  is  used  in  the  recursive  manner 
described  previously. 

7.3  Error  Covariance  of  Rayleigh  Estimator 

The  estimate  given  by  eq.  7.33  provides  a minimum 
mean-square  error  reconstruction  of  a quantized  Rayleigh 
variable.  The  performance  of  this  estimate  can  be 
determined  by  a calculation  of  the  estimation  covariance. 
A general  expression  for  this  covariance  was  derived  as 
eq.  3.10.  This  expression  can  be  solved,  for  a Rayleigh 
probability  density  function,  by  the  same  technique  that 
was  used  in  equations  5.36  to  5.37  for  a gaussian 
distribution.  Denoting  the  (i,i)th  component  of  the 
estimation  covariance  (the  only  non-zero  component  because 
of  the  recursive  nature  of  tno  estimation  procedure 
described  in  Sec.  7.2)  as  e. . , then 

li ' 
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or 


where 


e.  . =■ 
11 


1xi  pUJdx.^ 


j 1p(t)dxi 

J 


- 


(7.34) 


e.=„—  " (Vi+F"1)2 


11  e 


denom 


1 1 


(7.35) 


e =2 
num 


\fre  rf(z)  £ 
k = l 


IN+2 

LjU  g2k-2  (2k-1)  ! ! 
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k+1 


-exp{-z2 } £ gv  (k+1^ ! ! £ 


k-2 j+1 


k=0  k j = 0 2]  ( k-2 j+1)  1 ! 


(7.36) 


Ju 


e = r . . 

denom  n 


|]  92k(2k-i.n 

Viet£(z,  E ? 
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k-1 


-exp{-z2}  £ g(k-i) ! 1 £ 


k-2 j-1 


k = l 


j = 0 2-*  ( k-2  j-1 ) ! ! 


(7.37) 
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This  error  covariance  for  the  (i)th  component  is  related  to 
the  estimates  of  the  other  components,  tR  for  Mi,  by  mean 

by  eq . 7.27,  and  for  k=l ,/,... ,N , 


of  both  R. , as  given 
i 
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as  given  by  eq.  7.26.  The  expression  in  eq.  7.35  hence 
specifies  the  conditional  covariance  of  , conditioned  on 
the  other  estimates,  t^  . 

Because  the  components,  , are  random  variables,  the 
conditioning  can  be  removed  by  integrating  eq.  7.35  with 
respect  to  the  distribution  of  the  . Equivalently,  the 
dependence  on  the  t^  can  be  removed  by  integrating  with 
respect  to  and  , which  are  functions  of  the  t^  . The 
first  and  second  moments  of  and  can  be  calculated  as 


N 


E(V-r 


r.  . E { t . } 
13  3 


(7.38) 


E{ R7 } =-r . 

1 11 


N 

I 

3 = 1 

j^i 


r.  . C.  . 
i3  3i 


(7.39) 


_ N 

E{Hk}>2  .?  bfcj  k=l,2,...,N  (7.40) 

E,Hkls|1hkn|jCnjhkj  k=1'2 N (7-*D 

n^i  j^i 

It  can  be  shown,  by  means  of  lengthy  derivations,  that 


E{R.  } « E { H } 
l k 


(7.42) 
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and 


k 


} 


< 


E{R.2  } « } (7.43) 

for  k = 1 , 2 , . . . , N . Since  bot,  Ri  and  HR  consist  of  a sum  of 
many  random  variables,  it  can  be  argued  that  they  have  a 
gaussian  distribution.  Then,  under  this  assumption,  the 
two  moments  found  in  equations  7.38  and  7.39  are  sufficient 
to  characterize  this  distribution,  now  denoted  as  p{R^) . 
Hence,  the  unconditioned  error  covariance  can  be  found  from 

e = J eii  p(Ri)dP1  (7.44) 

However,  due  to  the  complicated  form  of  the  constituents  of 
this  expression,  the  error  covariance  cannot  be  evaluated 
for  general  cases.  Rather,  specific  cases  of  the 
estimation  covariance  must  be  considered  on  an  individual 
basis.  An  approximation  can  then  be  obtained  numerically. 

*/ . 4 Simulation  and  Restoration  of  PCM  C jded  Images 

To  determine  the  utility  and  effectiveness  of  the 

restoration  technique  developed  in  Sec.  7.2,  the  technique 

was  computer-simulated  and  then  applied  to  imaqes  which  had 

been  quantized  according  to  various  bit  assignments.  The 

i.'aqes  were  quantized  in  intensity  using  Max  decision  and 

restoration  levels.  These  levels  are  shown  in  table  7 1. 

They  were  calculated  by  means  of  eq.  2.3,  using  a 
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TABLE  7-1 


QUANTIZATION  INTERVALS  FOR  SIGNALS  WITH  A RAYLEIGH 
DISTRIBUTION  CHOSEN  ACCORDING  TO  A MINIMUM 
MEAN-SQUARE  ERROR  CRITERION 


Bits 

Quantization  Interval 

Restoration  Point 

0 

0.0  oo 

1.253 

1 

0.0  1.375 

0.829 

1.375  00 

1.920 

2 

0.0  0.822 

0.529 

0.822  1.420 

1.114 

1.420  2.127 

1.725 

2.127  OO 

2.529 

3 

0.0  0.499 

0.329 

0.499  0.825 

0.670 

0.825  1.135 

0.980 

1.135  1.453 

1.290 

1 

1 *453  : .800 

1.617 

1.800  2.208 

1.984 

2.208  2.760 

2 . *33 

1 

2.760  © o 

3.086 
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one-dimensional  Rayleigh  probability  density  function. 
Each  pixel  was  quantized  individually,  i.e.,  independently 
of  its  neighboring  pixels.  The  results  of  this  quantizing 
scheme,  when  applied  to  the  "girl"  image,  are  shown  in 
figures  7-5a  and  7-5c  for  one  bit  per  pixel  and  two  bits 
per  pixel  quantizers,  respectively. 

Figures  7-5a  and  7-5c  were  a posteriori  restored  by 
utilizing  ea.  7.33  with  N=5.  Using  this  value  for  N means 
that  the  restoration  of  one  pixel  is  achieved  by  utilizing 
the  information  that  is  available  about  four  other  pixels. 

A five-dimensional  correlated  Rayleigh  probability  density 

I 

function  is  thus  used  to  model  the  image  for  this 
restoration  procedure.  The  five  pixels  were  chosen  and 
ordered  as  shown  in  fig.  7-6.  These  particular  pixels 
allow  maximum  use  to  be  made  of  the  correlation  between 
adjacent  pixels.  Tne  5x5  correlation  matrix  of  these 
pixels  is  also  shown  in  fig.  7-6.  This  correlation  matrix 
is  a necessary  input  to  the  restoration  equation. 

The  restoration  provides  an  improvement  because  it 
permits  many  possible  output  intensity  levels  to  occur, 
even  though,  as  in  the  one  bit  case,  there  are  only  two 
input  levels  <*>.  The  restored  images  for  one  bit  and  for 
two  bits  are  shown  in  figures  7 — 5b  and  7-5d,  respectively. 
The  restoration  of  the  one  bit  quantized  image  results  in  a 
decrease  of  20.8%  in  mean-square  error.  The  restored  two 
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c.  Two  bits/pixel  PCM  d.  Restored  two  bits/pixel 


Figure  7-5.  Minimum  mean-square  error  reconstructions  of 
PCM  coded  images. 
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b.  Correlation  matrix  of  the  ordered  pixels 


Figure  7-6.  Typical  ordering  of  five  PCM  coded  pixels  for 
restoration,  and  their  corresponding  correlation  matrix. 


bit  image  has  28.6%  lower  mean-square  error  than  ;ts 
quantized  version.  Each  restored  image  exhibits  a subtle, 
but  noticeable,  visual  improvement  over  its  corresponding 
quantized  image.  The  restored  versions  appear  subjectively 
to  be  more  "real"  due  to  the  extra  intensity  levels  that 
result  from  the  restoration.  It  should  be  emphasized  that 
these  improvements  were  obtained  a posteriori,  utilizing 
only  the  correlation  matrix  of  a typical  sampled  image, 
together  with  the  quantized  image  to  be  restored. 


<*>  For  this  one  bit  case,  since  each  of  five  quantized 
pixels  can  be  in  one  of  two  possible  states,  there  are  32 
possible  output  values  for  the  restored  pixel. 
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CHAPTER  8 


RESTORATION  OF  BINARY  SYMMETRIC  CHANNEL  ERRORS 
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The  previous  chapters  have  presented  and  analyzed 
techniques  for  restoring  the  output  of  a quantizer  so  that 
the  result  more  accurately  matches  the  quantizer's  input. 
The  restorations  are  based  essentially  upon  exact  knowledge 
of  the  quantizer  output.  A similar,  but  more  difficult 
problem  results  when  the  quantizer  output  is  not  known 
exactly.  This  could  occur,  for  example,  when  the  quantizer 
output  is  transmitted  over  a noisy  channel.  The  first 
section  in  this  chapter  explores  the  effect  of  channel 
errors  on  the  restorations  derived  previously.  The  next 
section  examines  a technique  that  statistically  compensates 
for  the  effect  of  channel  errors. 


8.1  Effects  of  Channel  Errors  on  Quantized  Signals 


In  this  analysis,  channel  errors  are  assumed  to  arise 
in  the  context  of  a binary  symmetric  channel  (BSC)  [1]. 
The  characteristics  of  this  type  of  channel  are  shown  in 
fig.  8-1.  The  channel  is  discrete  and  memoryless  and  can 
be  specified  by  a transition  probability  assignment  P(j|k), 


for  j , k=0 , 1 , as 
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1-p 


p 


P={ P ( j I k)  } = 


Since  the  channel  is  memoryless,  the  probability  of  an 


output  sequence  z=  ( z , z2  , . . . , z^j ) , given  an  input  sequence 
x-(x  ,x  ,...,j^),  is  given  by 


1 ( Z | X ) = TT  P ( z . |x.  ) 

1=1  1 1 


(8.2) 


Based  on  this  definition,  a BSC  was  simulated  bv  means 


of  a computer,  with  the  channel  error  probability,  p, 


chosen  to  be  0.01.  The  simulated  channel  was  then  applied 


to  transform  coded  images.  The  three  "original"  images 


shown  in  fig.  5-7  were  zonal  transform  coded  in  16x16 


blocks,  as  described  in  Sec.  5.4.  The  quantized  transforr 


domain  components  were  encoded  by  assigning  each  a binary 


code  word.  The  resulting  sequence  of  binary  digits  was 


operated  on  by  the  simulated  channel.  The  error-corrupted 
bit  stream  was  then  either  decoded  directly,  as  shown  in 


figures  8-2a,  8-2c,  and  8-2e,  or  restored  by  the  techniaues 


of  Chapter  5 to  reduce  the  effects  of  the  quantizati< 


process  (see  fig.  8-3  for  a schematic  of  this  procedure). 


The  decoded  images  with  the  quantization  effects  reduced 


are  shown  in  figures  8-2b,  8-2d,  and  8-2f. 


Max  quantized 


Max  quantized 


e.  Max  quantized 


Figure  8-2.  Minimum  mean-square  error  restoration  of  Haar 
transformed  0.5  bit  zonal  quantized  images  transmitted 
through  a BSC  with  an  error  probability  = 0.01. 
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Bit  errors  in  transform  coding  that  arise  due  to  a 
binary  symmetric  channel  are  seen  to  result  in  an  emphasis 
of  the  block  structure  and  a subjective  error  that  extends 
over  the  entire  block.  This  latter  effect  occurs  because 
inverse  transforming  a block  containing  an  error 
distributes  this  error  over  all  the  resultant  image  domain 
components.  However,  what  is  important  to  note  from 
fig.  8-2  is  that  channel  errors  affect  quantized  and 
quantization-restored  images  identically.  The 
reconstruction  techniques  derived  in  Chapters  5,  6,  and  7 
are  thus  insensitive  to  channel  errors.  Since  they  provide 
visual  and  mean-square  error  improvements  in  noise-free 
cases,  they  can  be  utilized  equally  well  in  noisy 
environments . 

8.2  Reconstruction  of  Quantized  and  Transmitted  Signals 

The  previous  section  demonstrated  that  channel  errors 
do  not  adversely  affect  the  performance  of  the  restoration 
techniques  derived  earlier.  However,  these  techniques  do 
nothing  to  ameliorate  the  effects  of  the  channel  errors. 
This  is  because  the  fundamental  restoration  equation 
presented  in  eq.  3.6  was  derived  without  ar.y  consideration 
of  channel  structure.  By  including  the  channel  structure 
in  the  derivation,  the  resultant  restoration  technique  can 
simultaneously  reduce  the  effects  of  the  quantization 
process  and  mitigate  the  effects  of  channel  errors. 
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The  output  of  a data  source  (this  output  could  consist 
of  DPCM  samples,  PCM  samples,  or  transform  domain  samples) 
is  denoted  by  x=  ( ,x2  , . . . , XN ) and  described  by  a 
probability'  density  function  p(x).  The  reconstruction  of 
x,  after  x has  been  quantized  to  one  of  M regions  and 
channel-error  corrupted,  is  denoted  by  z^=(z^,z 
for  k=l,2,...,M  (refer  to  fig.  8-^)  . The  mean-square  error 
that  results  from  this  process  is 


2 Vk 


M M / 

$=Y  ^ P(mlk)  / (x-z  ) (x-z,  )Tp(x)dx  (8.3) 

k=l  m = l J R 

m 


This  error  can  be  minimized  by  proper  choice  of  the 
restoration  points,  z^ . Setting  the  partial  derivatives  of 
this  error  with  respect  to  z equal  to  zero  yields 


M 

Y P(ml  k) 

m=l 


Y P <m  | k) 

m=l 


xp(x) dx 


p (x) dx 


(8.4) 


for  k=l,2,...,M.  This  expression  is  the  noisy  channel 
version  of  eq.  3.6  and  provides  a minimum  mean-square  error 
estimate  of  the  input  to  a quantizer  based  on  the  output  of 
a noisy  channel,  the  characteristics  of  the  quantizer,  and 
the  a priori  statistics  of  the  input.  This  equation  is 
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also  a multidimensional  version  of  a result  first  derived 
in  [2].  For  a noiseless  channel,  the  channel  matrix  P 
becomes  the  identity  matrix  and  ea.  8.4  reduces  to  eq.  3.6. 
When  the  probability  volume  integrals  in  the  denominator  of 
ea.  8.4  are  all  equal,  which  is  approximately  true  for  Max 
quantization,  the  restoration  equation  simplifies  to 


M 


z = £ P(m | k) 


m=l 


I xp(x)dx 

Jr 

m 

/ p(x)dx 

J R 

m 


(3.5) 


or 


M 

ik=E  POnMOXjn 

m=l 

where  y is  given  by  eq.  3.6.  This  result  holds  for 

m 

maximum  output  entropy  quantizers  and  two-level  symmetrical 
quantizers,  and  is  approximately  correct  for  many  other 
types. 

A signal  that  has  been  quantized  and  then  transmitted 
over  a noisy  channel  can  thus  be  optimally  restored  by 
utilizing  eq.  8.4.  The  restoration  solutions  found  earlier 
for  gaussian,  laplacian,  and  Rayleigh  probability  density 
functions  (see  equations  5.23,  6.7,  and  7.33,  respectively) 
can  be  substituted  directly  into  eq.  8.4  once  the 

transition  matrix  for  the  channel  has  been  determined.  The 


resultant  estimator  can  then  be  used  to  restore  the  outputs 
of  transform,  DPCM,  and  PCM  coders  that  have  been  degraded 
by  channel  errors. 
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CHAPTER  9 


CONCLUSIONS,  AND  TOPICS  FOR  FUTURE  RESEARCH 

This  dissertation  has  described  a means  for 
reconstructing  quantized  signals  according  to  a minimum 
mean-square  error  criterion.  The  method,  which  is  entirely 
a posteriori,  is  based  on  a priori  statistical  information 
about  the  original,  unquantized  signals.  In  a broad  sense, 
the  reconstruction  technique  described  is  applicable  to  any 
coding  system  for  which  the  necessary  a priori  statistical 
information  is  available,  since  any  coding  system  can  be 
considered  as  a special  case  of  a quantizer.  In  other 
words,  assigning  one  of  a finite  number  of  code  words  to  a 
random  variable  is  equivalent  to  quantizing  that  variable 
to  one  of  a finite  number  of  intervals.  Further,  block 
encoding  a string  of  random  variables  is  the  sane  as 
quantizing  a random  vector  to  a generalized  region  in 
space . 

No  attempt  has  been  made  in  this  work  to  optimize  the 
location  or  choice  of  these  generalized  quantization 
regions,  but  only  to  utilize  arbitrary,  given  regions  in  an 
optimal  manner  to  obtain  a signal  reconstruction.  Finding 
an  optimum  set  of  regions  to  vector-quantize  a string  of 
random  variables  remains  an  unexplored  area,  but  one  in 
which  fruitful  research  can  be  conducted. 


195 


all  cases, 
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The  techniques  described  herein,  in 


achieved  a reduction  in  mean-square  quantization  error  when 


compared  to  a Max  restoration  scheme.  The  reason  for  this 
improvement  over  a method  which  purportedly  provides  the 
"minimum  mean-square  error,"  is  that  Max  s scheme  is 


memoryless  while  the  restoration  described  herein  requires 
memory.  The  memory  requirement  arises  in  the  form  of 
knowledge  of  the  quantized  samples  that  surround  a sample 
being  restored.  This  added  information , which  must  be 


stored  in  memory,  permits  the  improvement  in  mean-square 


error.  The  restoration  techniques  derived  in  this 


dissertation,  therefore,  represent  a generalization  of 


Max's  restoration  results  to  many  variables. 


From  an  information  theory  standpoint,  this  use  of 


memory  allows  a coding  system  to  operate  closer  to  the  rate 


distortion  theory  bound.  Also  utilizing  memory  to  perform 


source-encoding — equivalent  to  vector  quantization,  as 


discussed  previously — can  provide  an  even  closer  approach 


to  this  bound.  However,  in  this  work  a memoryless  encoding 


scheme  has  been  employed  with  its  implicit  assumption  of 


rectangular  quantization  regions. 


The  reconstruction  technique  developed  herein  has 


specifically  been  applied  to  three  common  types  of  image 


coding  systems — transform  coders,  DPCM  coders,  and  PCM 


coders.  These  coding  systems  have-  been  statistically 
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characterized  by  gaussian,  laplacian,  and  Rayleigh 
probability  density  functions,  respectively.  However,  to 
completely  characterize  the  coders  requires  correlated 
multidimensional  versions  of  these  density  functions. 
Heretofore,  multidimensional  versions  have  not  existed  for 
either  laplacian  or  Rayleigh  distributions.  This 
deficiency  motivated  the  general  technique  derived  in 
Chapter  4 for  generating  correlated  multidimensional 
density  functions  from  desired  marginal  distributions  and 
correlation  functions.  In  essence,  this  allows  the 
correlation  which  exists  between  the  variables  in  many 
coding  systems  to  be  utilized  in  a statistical  restoration 
of  their  coded  outputs. 

A determination  of  the  minimum  mean-square  error 
restoration  point  for  a quantized  vector  has  often  proven 
difficult,  because  a multidimensional  integration  of  a 
complicated  probability  function  is  required.  However, 
this  difficulty  has  been  surmounted  by  a novel,  recursive 
approach.  This  approach  permits  a reduction  to  only  a 
single  integration  which  can  then  be  evaluated 
analytically.  The  solutions  which  resulted  from  this 


approach  have 

been 

applied 

to  quantized 

images 

, and 

only 

one  recursion 

of 

the  technique  has 

been 

found 

to  be 

necessary  for 

an 

image 

restoration . 

The 

restoration 

procedure  can  thus  be  readily  implemented. 
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A decrease  in  mean-square  error  was  obtained  in  all 


I 
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simulations  of  this  restoration  procedure.  However,  a 
corresponding  subjective  improvement  was  not  always 
obs^.ved.  This  ancillary  result  substantiates  an 
observation  made  by  many  others:  lower  mean-square  error 
does  not  always  correspond  to  a visual  improvement, 
repeating  the  image  restorations  in  the  context  of  a 
weighted  error  criterion,  with  the  weighting  chosen 
according  to  characteristics  of  the  human  visual  system, 
resulted  in  both  an  analytical  and  a subjective 
improvement.  It  is  suggested  that  further  image  coding 
reconstructions  be  performed  with  respect  to  a weighted 
error  criterion. 

A fundamental  limitation  imposed  on  the  images 
restored  in  this  dissertation  was  an  assumption  of 
statistical  stationar ity . This  restriction  is  very  basic, 
because  the  reconstruction  techniques  herein  are  completely 
dependent  on  the  choice  of  a statistical  model  for  the 
underlying  random  process.  Since  images  have  been  found  to 
be  inherently  nonstationary,  an  assumption  of  stationarity 
limits  the  reconstruction  performance.  An  adaptive 
technique  could  remove  the  limitation  and  improve  the 
performance.  This  has  not  yec  been  investigated,  but  could 
prove  to  be  a productive  area  of  research. 

Another  problem  that  was  discussed  only  briefly  in 
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this  work  concerns  the  reconstruction  of  quantized  signals 
in  the  presence  of  channel  noise.  This  is  a much  more 
difficult  problem  than  noise-free  quantization  restoration 
and  warrants  further  study. 


